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Abstract 

This  study  is  based  on  the  merger  of  two  separate  theories  to  further  the  efficiency 
with  which  joined-wing  structural  models  are  designed.  The  first  theory  is  Geometrically 
Exact  Beam  Theory  (GEBT).  GEBT  is  a  small  strain  beam  theory  which  is  capable 
of  accurately  capturing  the  geometric  bend-twist  coupling  in  beam  elements  that  are 
experiencing  large  global  deformations.  This  is  crucial  to  the  joined-wing  problem  as  it 
is  geometrically  nonlinear.  The  second  theory  concerns  Equivalent  Static  Loads  (ESL). 
These  ESL  consist  of  a  load  vector  that  produces  the  same  nodal  displacements  and 
rotations  as  those  computed  from  a  pure  nonlinear  analysis.  The  ESL  displacements 
and  rotations  are  then  used  to  calculate  ESL  stresses.  By  merging  these  two  theories 
into  a  single  structural  optimization  effort,  computational  cost  is  reduced  by  orders 
of  magnitude  when  compared  to  purely  nonlinear  response  optimization  efforts.  It  is 
shown  that  the  final  design  obtained  by  the  optimization  is  the  same  for  both  types  of 
analysis.  The  final  result  is  a  much  simpler  model  than  a  detailed  finite  element  model 
of  the  joined-wing  aircraft  that  can  be  optimized  without  significant  loss  in  fidelity  in 
a  fraction  of  the  time  required  for  a  single  nonlinear  response  optimization  cycle  using 
finite  element  analysis. 
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Structural  Optimization  of  Joined-Wing  Beam  Model 
with  Bend-Twist  Coupling 
Using  Equivalent  Static  Loads 

I.  Problem  Statement 

The  joined-wing  configuration  presented  by  Wolkovich  in  his  1986  paper  has  at¬ 
tracted  much  interest  due  to  the  advantages  presented  by  such  a  configuration. 
Chief  among  these  advantages  are  an  increase  in  stiffness  when  compared  to  a  classic 
wing  and  tail  arrangement  for  a  given  weight  when  designed  correctly  [3].  Optimization 
efforts  on  joined-wing  aircraft  have  steadily  increased  in  complexity  from  fully  stressed 
design  [4]  to  efforts  that  incorporate  the  linear  [5]  and  nonlinear  [6]  dynamic  response 
using  equivalent  static  loads.  All  of  the  aforementioned  work  has  been  done  using  finite 
element  analysis  as,  until  recently,  beam  theory  was  unable  to  fully  capture  the  nonlinear 
geometric  bend-twist  coupling  observed  in  the  joined  wing  [7-10]. 


Figure  1.1:  Boeing  Sensorcraft  design  as  proposed  to  Air  Force  Research  Laboratory 
(AFRL) 

The  motivation  for  this  research  effort  stems  directly  from  the  immense  computa¬ 
tional  effort  required  to  solve  nonlinear  problems  with  thousands  of  degrees  of  freedom. 
An  example  of  a  complex  joined-wing  model  that  to  date  has  only  been  analyzed  using 
finite  element  analysis  is  the  Boeing  Sensorcraft,  shown  in  Fig.  1.1.  This  design  has 
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been  proposed  to  the  Air  Force  Research  Laboratory  (AFRL)  as  a  possible  solution  to 
a  high  altitude,  long  range,  long  endurance  surveillance  platform  which  will  leverage  the 
unique  geometry  of  the  joined-wing  design  in  the  deployment  of  its  sensor  arrays.  Such 
a  complex  design  would  almost  certainly  benefit  from  an  optimization  effort  which  could 
include  sizing  of  the  wing  box  dimensions,  the  angles  in  both  planform  and  head  on 
views  of  the  joined-wing,  aerodynamic  loads,  and  sensor  effectiveness  based  on  displace¬ 
ment  of  the  wing  from  a  nominal  configuration  in  flight.  Based  on  the  cost  of  a  single 
static  nonlinear  analysis  on  a  model  of  the  Boeing  design,  it  is  clear  that  an  optimization 
that  considers  all  of  these  design  variables  in  a  nonlinear  response  optimization  would 
be  unwieldy  at  best  and  intractable  at  worst.  In  an  effort  to  reduce  this  computational 
overhead,  this  research  explores  the  use  of  nonlinear  structural  analysis  on  a  beam  model 
of  the  joined-wing  in  a  linear  optimization  routine  using  equivalent  static  loads  in  lieu  of 
the  nonlinear  static  response  due  to  tip  loads. 

Early  studies  on  the  joined-wing  show  that  bending  due  to  a  vertical  tip  load  (with 
respect  to  the  longitudinal  axis  of  the  vehicle)  does  not  occur  in  the  plane  normal  to  the 
wing  chord,  but  rather  in  the  direction  perpendicular  to  the  plane  created  by  the  axes 
of  the  fore  and  aft  wings  [3].  This  out-of-plane  bending  is  then  coupled  with  a  twisting 
of  the  wing  structure  due  to  the  geometry  of  the  wings. 

Simpler,  classic  wing  box  designs  were  often  modeled  as  an  Euler-Bernoulli  beam 
and  then  sized  using  fully  stressed  design  techniques  for  isotropic  materials.  As  comput¬ 
ers  became  more  powerful  in  the  90’s,  finite  element  analysis  allowed  for  accurate  cal¬ 
culations  involving  complex  cross-sections  with  non-isotropic  material  properties.  While 
the  results  from  finite  element  analysis  are  valuable  in  the  design  of  complex  structures, 
they  are  still  time  consuming.  A  simple  yet  complete  beam  theory  would  drastically 
reduce  computational  effort  via  a  reduction  in  the  number  of  elements,  but  an  adequate 
beam  model  did  not  exist  before  Hodges’  work  on  his  intrinsic  beam  theory  [11],  With 
continuing  refinement  and  validation,  a  number  of  structural  specialists  have  developed 
successful  FEM  beam  formulations  based  on  Hodges’  theory  with  an  eye  towards  heli¬ 
copter  rotor  design  applications  [12-19].  Recent  strain-based  methods  have  been  applied 
to  the  analysis  and  design  of  various  equivalent-beam  wing  structures  [20,  21],  Geo- 
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Figure  1.2:  Planform  and  view  looking  forward  of  joined- wing  consisting  of  eight  beam 
elements. 

metrically  Exact  Beam  Theory  (GEBT)  is  an  alternate  beam  element  formulation  with 
mixed  displacement  and  force  variables  (as  opposed  to  strain  variables).  Hodges’  in¬ 
trinsic  beam  theory  incorporates  all  possible  geometric  nonlinear  coupling  combinations 
including  geometric  bend-twist  coupling  reported  by  Blair  and  Stritz  [10].  That  devel¬ 
opment  motivated  Yu  in  the  broader  development  of  GEBT  with  the  goal  of  supporting 
gradient-based  design  optimization  of  slender  aeroelastic  wing  concepts. 

At  the  same  time,  Park  et  al.  developed  an  optimization  scheme  for  linear  [5] 
and  nonlinear  [6]  joined-wing  structures  based  on  equivalent  static  loads.  Simply  put, 
the  displacements  from  the  nonlinear  analysis  are  used  to  develop  a  load  vector  that 
yields  the  same  set  of  displacements  when  applied  to  a  linear  stiffness  matrix.  This 
two-step  method  starts  with  the  converged  solution  of  the  nonlinear  structural  system, 
then  applies  a  linear  approximation  technique  to  the  current  design.  The  equivalent  load 
vector,  or  equivalent  static  loads  (ESL),  are  then  used  in  the  linear  response  optimization 
until  convergence  is  achieved.  The  entire  process  then  repeats  until  the  equivalent  static 
loads  for  two  consecutive  nonlinear  evaluations  have  converged. 

The  effort  reported  in  this  thesis  document  is  a  combination  of  GEBT  and  ESL  as 
developed  by  Yu  and  Park  respectively.  A  current  limitation  of  this  implementation  of 
GEBT  is  the  inability  to  accommodate  loads  applied  to  the  wing  tip  beyond  the  wing 
joint.  Therefore,  only  a  design  where  the  wings  are  joined  at  the  wing  tip  is  considered. 


(b)  View  looking  forward 
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Load  conditions  will  include  only  static  vertical  loads  applied  to  the  joint  at  the  wing 
tip. 

The  goals  of  this  research  are  first  to  develop  a  pilot  application  of  the  merger  of 
GEBT  and  ESL.  This  merger  leads  to  the  second  goal,  which  is  to  demonstrate  the  util¬ 
ity  of  GEBT  in  a  computationally  efficient  design  optimization  process  of  a  joined-wing 
concept  using  the  method  of  equivalent  static  loads.  The  simplest  possible  realization 
which  will  exhibit  geometric  bend-twist  coupling  in  a  joined  wing  is  shown  in  Fig.  1.2. 
Second,  show  that  a  low  degree  of  freedom  model  of  the  joined  wing  is  capable  of  cap¬ 
turing  geometric  bend/twist  coupling  exhibited  in  the  joined-wing  design.  This  thesis 
reports  on  the  success  of  this  demonstration.  GEBT  will  continue  to  be  developed  in  a 
series  of  steps  that  increase  its  usefulness  in  a  computationally  efficient  aeroelastic  design 
optimization  process  applied  to  highly  flexible  concepts. 
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II.  Background 


THIS  research  relies  on  three  main  bodies  of  work  as  foundational  material.  The 
first  is  the  joined-wing  design  itself.  The  second  is  Geometrically  Exact  Beam 
Theory  and  its  associated  Variational  Asymptotic  Method  and  Variational  Asymptotic 
Beam  Section  theories  for  analysis  of  beam  elements.  Third,  Equivalent  Static  Loads  will 
be  discussed  as  a  means  of  reducing  computational  effort  to  solve  nonlinear  structural 
optimization  problems. 

2. 1  The  Joined-  Wing  Aircraft 

The  joined-wing  aircraft  design  is  not  a  new  concept  within  the  history  of  manned 
flight.  The  novelty  of  the  joined-wing  is  that  the  tools  necessary  to  analyze  the  complex 
structural  and  aerodynamic  properties  were  not  available  until  the  mid  to  late  1980’s. 
History’s  first  known  successful  flight  of  a  joined-wing  design  was  made  by  Reinhold 
Platz  in  1922  [1],  The  Platz  joined-wing  glider  featured  two  wings  joined  at  the  tip,  with 
the  split  upper  forward  wing  being  used  to  control  both  pitch  and  roll  by  the  pilot.  The 
lower  aft  wing  was  fixed  and  supplied  most  of  the  lift  as  shown  in  Fig.  2.1. 

Other  designers  experimented  with  the  joined-wing,  but  it  wasn’t  until  Wolkovitch 
was  granted  his  patent  in  1982  and  presented  an  overview  of  the  design’s  merits  that 
aircraft  designers  began  to  consider  the  configuration  as  feasible  [2],  Wolkovitch’s  patent 
application  drawing  is  shown  in  Fig.  2.2.  Several  advantages  are  attributed  to  the  joined- 
wing  design  when  compared  to  a  conventional  aircraft  design,  including  light  weight,  high 
stiffness,  low  induced  and  parasite  drag,  direct  lift  control,  and  good  stability  and  control 
[3], 

These  advantages  have  since  prompted  much  research  into  the  joined-wing  design. 
Structural  optimization  of  the  joined-wing  design  was  investigated  by  Gallman  and  Kroo 
in  1996.  They  used  two  structural  design  methods:  one  a  fully  stressed  design  and  the 
second  a  minimum  weight  design  [4].  They  concluded  that  the  fully  stressed  design  us¬ 
ing  nonlinear  analysis  tools  approximates  the  minimum  weight  structure  very  well  with  a 
time  savings  in  computational  effort.  More  recently,  the  Air  Force  Research  Laboratory 
(AFRL)  has  shown  interest  in  the  joined-wing  design  concept  as  a  surveillance  platform 
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Figure  2.1:  1922  Platz  glider  [1] 


with  long  range  and  endurance  as  key  performance  parameters.  Blair  and  Canfield  devel¬ 
oped  an  integrated  design  method  for  the  joined-wing  aimed  at  merging  the  aerodynamic 
and  structural  designs  into  one  all  encompassing  process  [22],  They  concluded  that  the 
geometric  nonlinearities  inherent  in  the  joined-wing  require  a  nonlinear  analysis  for  loads 
that  result  in  large  deformations.  In  that  vein,  numerous  studies  have  been  performed  on 
the  joined-wing  using  the  finite  element  method,  with  models  having  as  many  as  3,000 
elements  being  considered  for  a  single  nonlinear  analysis  on  a  half  span  model  [6]. 

2.2  Geometrically  Exact  Beam  Theory  (GEBT) 

The  most  fundamental  aspect  of  elasticity  is  Hooke’s  Law,  which  in  ID  form  is 

F  =  —kx  (2.1) 

where  k  is  the  stiffness  of  the  element,  x  is  the  displacement  of  the  element,  and  F  is 
the  axial  force  applied  to  the  element.  When  expanded  beyond  the  ID  axial  problem, 
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(a)  Isometric  view 


(b)  Planform  view 


(c)  Side  view 


Figure  2.2:  Joined-wing  aircraft  design  proposed  by  Wolkovitch  in  his  patent  applica¬ 
tion.  [2] 
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application  of  Euler-Bernoulli  beam  theory  considers  only  a  4x4  stiffness  matrix  due  to 
the  selection  of  displacement  and  slope  at  the  two  nodes  as  the  degrees  of  freedom  un¬ 
der  consideration.  The  assumption  that  planar  surfaces  remain  in  plane  and  normal  to 
the  longitudinal  axis  under  deformation  allows  shear  deformation  to  be  neglected.  Timo¬ 
shenko  expands  on  the  Euler-Bernoulli  beam  theory  by  considering  planar  surfaces  which 
remain  plane,  but  not  normal  to  the  longitudinal  axis  of  the  beam  element.  Inclusion  of 
the  angle  of  shearing  in  the  two  directions  of  the  beam  element  cross  section  increases 
the  degrees  of  freedom  to  6,  resulting  in  a  6x6  stiffness  matrix. 

With  an  eye  towards  helicopter  rotor  blade  design,  Hodges  began  an  investigation 
into  beam  theory  in  the  mid  1980’s  that  would  eventually  result  in  a  generalized  beam 
theory  that  is  not  burdened  by  several  assumptions  that  are  often  required  when  con¬ 
ducting  analysis  of  long,  slender  beams.  He  presented  his  mixed  variational  formulation 
beam  theory  in  1989  and  notes  that  his  work  is  the  assemblage  of  three  sets  of  previous 
works  [11]. 

First  is  a  kinematical  basis  for  deformed  beams.  Danielson  and  Hodges  used  the 
polar  decomposition  theorem  to  obtain  an  accurate  yet  compact  expression  for  the  strain 
in  a  beam  or  rod  experiencing  large  deflections  [23].  The  polar  decomposition  theorem 
facilitates  this  by  decomposing  the  total  deformation  of  any  element  into  a  pure  strain 
and  a  pure  rotation.  Danielson  and  Hodges  concluded  that  strain  components  are  explicit 
functions  of  position  within  the  cross  section  and  implicit  functions  of  the  length  along 
the  beam  axis  X\.  There  are  a  total  of  seven  implicit  functions,  dependent  on  strain 
(711)  and  transverse  shear  (712  and  713)  of  the  reference  line,  bending  strains  («2  and 
K3),  and  warp  amplitude  (i/j).  In  1988,  Danielson  and  Hodges  developed  both  an  intrinsic 
and  explicit  version  of  their  beam  theory  [24],  The  intrinsic  beam  theory  is  for  static, 
nonlinear  beams  and  is  expressed  in  terms  of  the  generalized  strains  described  in  Ref.  [23] 
and  is  applicable  when  generalized  forces  and  moments  applied  at  the  ends  of  the  beam 
are  known.  The  explicit  theory  by  contrast  is  applicable  when  the  generalized  transverse 
and  rotational  displacements  are  known  in  conjunction  with  the  warp  amplitude  (iq,  0U 
and  if)  respectively).  Thus,  Danielson  and  Hodges  put  forward  theories  that  addressed 
the  global  kinematics  of  deformed  beams  using  Cartesian  tensors  and  matrix  notation, 


eliminating  the  need  for  curvilinear  coordinate  systems  and  reducing  the  size  of  the 
equations  required. 

Second,  Atilgan  and  Hodges  expanded  on  Parker’s  1979  work  in  Ref.  [25].  Parker 
showed  that,  using  a  perturbation  procedure  based  on  scaling  of  the  axial  variations,  St. 
Venant’s  semi-inverse  solutions  for  bending  and  torsion  of  rods  and  beams  could  be  used 
to  describe  beams  undergoing  large  displacements  and  rotations.  Parker  also  stated  that 
the  classical  solution  for  beam  flexure  is  composed  of  two  parts:  a  bending  solution  and 
a  varying  curvature  along  the  axis  of  the  beam.  The  use  of  asymptotic  analysis,  such  as 
that  used  by  Parker,  allowed  Atilgan  and  Hodges  to  conclude  a  linear  2D  cross-section 
deformation  analysis  can  be  used  to  determine  the  appropriate  stiffness  for  use  in  the 
ID  nonlinear  beam  global  analysis  [26]. 

Third,  Reissner  derived  strain  measures  from  internal  beam  forces  and  virtual  work 
using  an  intrinsic  analysis  [27,  28].  This  intrinsic  analysis  supports  the  later  work  by 
Danielson  and  Hodges  which  was  based  strictly  on  kinematics  [23].  The  agreement  be¬ 
tween  the  two  approaches  indicates  that  the  general  equations  to  describe  strain  measures 
are  now  available  to  the  engineering  community. 

Building  upon  the  three  principles  just  described,  Hodges  applies  a  variational 
approach  using  Hamilton’s  principle  that  concludes  with  a  mixed  variational  formulation 
for  a  geometrically  exact  beam  theory  [11],  The  advantages  with  this  formulation  lie  in 
the  compact  matrix  form  in  which  the  theory  is  expressed,  along  with  the  absence  of  any 
approximations  in  the  geometry  of  the  deformed  elements  along  the  reference  line  (i.e. 
beam  axis)  or  the  cross  section.  Hodges  also  notes  that  the  mixed  formulation  allows 
the  use  of  simple  shape  functions  consisting  of  constant  values  for  held  variables  along 
the  beam  element.  Discontinuities  in  held  variables  between  elements  are  also  possible, 
which  eliminates  the  requirement  of  numerical  quadrature  over  the  elements.  Instead, 
a  lxl  quadrature  is  performed  at  the  center  of  each  element  to  obtain  an  approximate 
stiffness,  mass,  etc. 

Hodges’  theory  has  been  expanded  to  use  the  Timoshenko  approach  by  including 
the  shear  terms,  but  still  relies  on  the  presence  of  small  strains  while  not  restricting  dis¬ 
placements  and  rotation  variables.  This  trait  is  due  to  the  derivation  of  the  theory  from 
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geometrically  nonlinear,  three  dimensional  elasticity  theory.  In  the  last  two  decades  the 
geometrically  exact  canonical  equations  of  motion  for  beams  have  been  presented  sep¬ 
arately  by  Borri  and  Mategazza,  Hodges,  and  Bauchau  and  Kang  [11,  29,  30].  This 
body  of  work  lends  itself  to  expressing  energy  in  terms  of  variables  compatible  with  large 
global  displacements  and  rotations.  Hodges’  Geometrically  Exact  Beam  Theory  (GEBT) 
requires  the  cross-sectional  clastic  constants,  El,  GJ,  etc.  as  the  inputs  to  the  equations 
used  for  solving  for  the  displacements  and  rotations.  While  this  approach  allows  for  the 
use  of  complex  beam  structures,  it  leaves  the  user  with  the  task  of  providing  the  cross- 
sectional  constants  before  the  calculation  can  proceed.  In  GEBT,  the  stiffness  matrix 
for  a  prismatic  isotropic  material  is 

^  EA  0  0  0  0 

0  GAXX  0  0  0 

0  0  GAZZ  0  0 

00  0  GJ  0 

0  0  0  0  EIXX 

y  0  0  0  0  0 

where  EA  is  the  axial  rigidity,  EIXX  is  the  bending  rigidity  about  the  global  X-axis, 
EIZZ  is  the  bending  rigidity  about  the  global  Z-axis,  GJ  is  the  torsional  rigidity,  GAXX  is 
the  shear  rigidity  in  the  global  X-direction,  and  GAZZ  is  the  shear  rigidity  in  the  global 
Z- direction.  The  ID  displacement  vector  e  is  defined  as 

e  =  |_7n  27i2  27i3  «i  k2  k3\t  (2-3) 

where  7n  is  the  extensional  strain  of  the  reference  line,  2712  and  2713  are  the  transverse 
shear  strain  measures,  and  /q  are  the  twist  and  bending  generalized  strain  measures  [12]. 
The  force  vector  applied  to  the  beam  is 

F  =  |_F,  F2  F3  Mi  M2  M3\t  (2.4) 


0 

0 

0 

0 

0 

El, 


\ 
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where  iq  is  the  force  in  the  X,  Y,  or  Z  direction,  respectively,  and  M*  is  the  moment 
about  the  X,  Y,  or  Z-axis.  Using  the  mixed  formulation  in  GEBT,  the  ID  constituitive 
relation  for  isotropic  materials  is 


(  \ 

F\ 

EA 

0 

0 
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0 
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7n 
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GAXX 
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0 
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«  3 

(2.5) 


There  are  at  present  two  implementations  of  GEBT  in  code.  Yu  has  made  avail¬ 
able  a  Fortran  90/95  based  version  that  as  of  this  writing  is  only  capable  of  addressing 
linear  static  problems.  The  GEBT  implementation  used  in  this  research  is  a  prototype 
Mathematica  code  that  allows  nonlinear  analysis  of  the  joined-wing.  The  inputs  to  the 
GEBT  implementation  are  the  applied  load,  the  number  of  elements,  the  total  length  of 
all  beam  elements,  and  the  cross  sectional  properties  for  each  element  described  in  Eq. 
(2.5).  The  GEBT  implementation  returns  forces,  moments,  displacements,  and  rotations 
for  each  element,  as  well  as  the  reaction  forces  at  the  constrained  nodes  (fore  and  aft 
wing  roots  in  the  case  of  this  research  effort).  The  forces,  moments,  displacements,  and 
rotations  are  reported  at  the  center  of  each  beam  element,  ft  is  possible  to  extrapo¬ 
late  nodal  displacements  and  rotations  by  starting  at  the  wing  root  and  using  a  linear 
interpolation 

U node  i  W node  i—1  Y  2 (Uelem  i  U node  i— l)  (2.6) 


where  unode  %  is  the  nodal  displacement,  unode  i- 1  is  the  nodal  displacement  at  the  previous 
node,  and  ueiem  %  is  the  displacement  at  the  mid  section  of  the  beam  element. 


2.3  Variational  Asymptotic  Beam  Sectional  Analysis  (VABS) 

The  Variational  Asymptotic  Beam  Sectional  Analysis  (VABS)  is  a  method  which 
uses  cross  sectional  properties  in  a  2D  analysis  based  on  the  variational  asymptotic 
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method  (VAM).  Some  of  the  background  on  VABS  is  described  in  Section  2.2.  The  ID 
nonlinear  beam  analysis  and  the  2D  cross-sectional  analysis  were  derived  in  the  context 
of  the  VAM  developed  by  Berdichevskii  in  1976  [31].  Berdichevskii  dealt  with  nonhomo- 
geneous  anisotropic  beams.  Hodges  et  al.  extended  the  work  to  include  prismatic  beams 
[7].  The  application  of  VAM  to  beam  theory  requires  that  the  3D  strain  energy  of  a 
beam  is  asymptotically  reduced  in  two  distinct  steps  to  yield  a  ID  beam  strain  energy 
equation.  This  ID  strain  energy  equation  is  decoupled  from  the  2D  cross-sectional  anal¬ 
ysis,  allowing  the  complex  3D  problem  to  be  split  into  two  separate  problems.  Popescu 
and  Hodges  offer  a  definition  of  ‘asymptotically  correct’: 

By  asymptotically  correct,  we  mean  that  an  expansion  of  the  approximate 
solution  in  terms  of  a  small  parameter  agrees  with  a  similar  expansion  of  the 
exact  solution  up  to  a  certain  order  in  the  small  parameter  [32], 

Yu  et  ah  expanded  Hodges’  work  to  handle  Timoshenko  like  modeling  of  initially  curved 
and  twisted  beams  to  address  short  wavelength  motions  caused  by  shear  effects  [9].  At 
the  same  time,  he  updated  an  engineering  software  package  known  as  VABS  III,  hereafter 
called  VABS,  to  include  the  latest  refinements  to  Hodges’  model.  VABS  is  a  computer 
program  based  on  the  VAM  which  allows  users  to  analyze  composite  beam  cross  sections 
accurately  and  quickly  using  finite  element  methods  with  the  processing  power  found  in 
the  average  desktop  computer  [33]. 

An  input  hie  is  required  that  describes  the  cross  section,  material  properties,  coordi¬ 
nate  transformation,  and  forces  and  moments  applied  to  the  element  under  consideration 
[34],  The  cross  section  is  defined  by  a  set  of  nodes  that  refer  only  to  points  within  the 
cross  section,  and  have  no  relation  at  all  to  the  nodes  depicted  in  Fig.  1.2.  For  this 
research,  VABS  is  used  to  supply  highly  reliable  calculations  of  the  nonlinear  stresses 
resulting  from  the  forces  and  moments  calculated  in  the  nonlinear  analysis  by  GEBT  [8] . 

2-4  Calculation  of  Equivalent  Static  Loads 

The  geometric  nonlinearity  of  the  joined-wing  leads  the  designer  down  a  difficult 
path.  A  linear  optimization  is  not  likely  to  be  adequate,  or  perhaps  even  safe  as  noted 
by  Kim  et  al.  [6].  During  cruise  and  maneuvering  flight,  Kim  et  al.  showed  that  the 
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maximum  tip  von  Mises  stress  in  a  geometric  nonlinear  analysis  could  be  as  much  as 
10  times  the  maximum  von  Mises  stress  in  a  linear  analysis,  while  displacements  in  the 
nonlinear  case  were  5  times  the  maximum  displacements  in  the  linear  analysis.  From 
this  it  is  clear  that  a  linear  analysis,  while  fairly  simple  to  implement,  is  not  sufficient 
when  considering  design  of  a  structure  that  exhibits  highly  nonlinear  responses. 

The  alternative  to  the  linear  response  optimization  has  been  a  full  nonlinear  re¬ 
sponse  optimization.  The  major  drawback  to  any  nonlinear  problem  lies  in  the  large 
computational  effort  involved  in  calculating  a  solution.  This  increase  in  effort  is  due 
to  the  numerical  approach  of  convergence  to  a  solution  using  multiple  intermediate  lin¬ 
earizations.  Sensitivity  calculations  in  the  optimization  contribute  to  the  total  effort  as 
well.  Analytic  sensitivity  calculation  and  the  finite  difference  method  are  the  two  com¬ 
mon  approaches  to  obtaining  sensitivities.  Although  computationally  expensive,  finite 
difference  techniques  are  easy  to  implement.  Analytical  evaluations  of  sensitivity  have 
advantages  and  drawbacks  that  are  opposite  of  the  finite  difference  method.  Analytical 
means  are  cheap  in  terms  of  computational  effort,  but  are  often  either  very  complex  or 
simply  undetermined  in  large,  complex  nonlinear  problems.  With  the  use  of  the  finite 
difference  method,  it  follows  that  any  increase  in  the  number  of  design  variables  will 
directly  affect  the  length  of  time  required  to  perform  each  iteration  in  the  optimiza¬ 
tion  process.  In  an  effort  to  circumvent  large  numbers  of  nonlinear  operations,  a  linear 
approximation  to  the  nonlinear  response  optimization  problem  has  been  proposed. 

Lee,  et  al.  define  equivalent  loads  (EL)  as  loads  for  linear  analysis  which  generate 
the  same  response  fields  as  those  of  nonlinear  analysis  [5].  In  the  optimization  problem  in 
this  research,  equivalent  static  loads  (ESL)  are  used  to  define  a  linear  approximation  of 
the  displacement  field  resulting  from  a  nonlinear  response  to  an  applied  load.  The  ESL 
are  then  passed  to  the  optimization  routine  allowing  a  linear  response  optimization  to 
be  performed,  resulting  in  less  computational  effort  required  to  determine  a  solution  to 
the  optimization  problem.  After  convergence  of  the  linear  approximation,  the  design  is 
reevaluated  using  a  nonlinear  analysis  on  the  linear  approximation  design  and  checked  for 
convergence.  When  the  equivalent  static  loads  are  no  longer  changing  from  one  iteration 
to  the  next,  the  optimization  problem  has  converged. 
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2-4-1  ESL  for  Displacement  Constraints.  The  first  step  in  calculating  an  equiv¬ 
alent  static  load  is  to  perform  a  nonlinear  analysis 

K(b)zn  =  f  (2.7) 

where  K(b )  is  the  nonlinear  stiffness  matrix,  zn  is  the  nonlinear  nodal  displacement 
vector,  and  f  is  the  applied  load  vector.  Due  to  the  dimensional  reduction  of  the  3D 
beam  problem  when  using  VAM,  the  cross  sectional  analysis  becomes  a  linear  2D  prob¬ 
lem  involving  the  cross  sectional  stiffness  properties  [34],  This  is  different  than  the 
approach  described  by  Lee  and  Kim  where  the  stiffness  matrix  is  a  function  of  both  the 
design  variables  b  and  the  displacement  vector  zn  [5,  6].  Equivalent  static  loads  ff  for 
displacements  are  computed  from 


fzeq  =  KL(b)zn  (2.8) 

where  K^ip)  is  the  linear  stiffness  matrix,  zn  is  from  the  nonlinear  analysis  from  Eq. 
(2.7). 


2-4-2  ESL  for  Stress  Constraints.  Similarly,  the  equation  for  the  equivalent 
static  load  (ESL)  for  stresses  f*  is 


fZ,  =  KL(b)z°n  (2.9) 

where  is  the  displacement  vector  calculated  from  the  nonlinear  stress  result  described 
[5].  In  the  research  presented  here,  an  alternate  means  of  calculating  the  equivalent 
stresses  is  used.  This  alternate  approach  uses  the  same  equivalent  static  loads  vector 
calculated  for  the  displacement  constraint  case  and  is  presented  in  MSC  Nastran’s  im¬ 
plementation  of  the  ESL  routine  [35].  The  equivalent  static  load  for  displacement  does 
not  produce  the  same  set  of  stresses  as  the  nonlinear  analysis,  so  the  stress  ratio  correc¬ 
tion  from  Lee,  et  a.1.  is  retained  [5].  The  stress  ratio  o  is  defined  as 

vnl-vm 

OL  —  - 

a  L  —  VM 


(2.10) 


14 


which  is  calculated  prior  to  entering  the  linear  response  optimization.  The  von  Mises 
stress  from  nonlinear  analysis  is  ctnl-vm^  and  the  von  Mises  stress  from  the  linear 
approximation  using  equivalent  static  loads  is  ctl-vm-  In  the  optimization  routine,  it 
is  applied  to  the  linear  stress  calculation  using  linear  Euler-Bernoulli  beam  theory.  The 
stress  vector  approximation  cr L_VM  is  defined  as 

l—vm  =  a  ^l—vm  (2-11) 

where  cr L  is  the  stress  vector  calculated  using  as  the  load  before  the  correction 
factor  ol  is  applied.  The  vector  cL_VM  is  used  to  determine  the  von  Mises  stress  in 
each  element.  These  values  are  then  used  in  stress  constraint  violations  such  that 

aL-VM  S 

where  crVM-allowable  js  the  von  Mises  stress  constraint  defined  in  the  optimization 
problem. 
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III.  Methodology 


3.1  Model  Definition 

The  model  is  described  in  Cartesian  coordinates  using  a  global  right  hand  coordi¬ 
nate  system  with  the  wing  tip  as  the  origin.  The  stream  wise  direction  (the  direction 
the  airflow  would  go  if  the  vehicle  were  in  flight)  is  along  the  positive  global  X-axis. 
The  positive  Y-axis  is  directed  out  the  span  through  the  wing  tip,  and  is  normal  to  the 
longitudinal  axis  of  the  aircraft.  The  positive  Z-axis  is  directed  in  the  upward  direction. 
The  beam  elements  are  oriented  so  that  the  area  exposed  to  the  free  stream  flow  is  the 
height  of  each  beam  multiplied  by  its  length  as  shown  in  Fig.  3.1. 


Figure  3.1:  View  looking  aft  showing  beam  element  orientation  (not  to  scale) 

Wolkovitch  states  that  the  effect  of  joint  location  has  a  great  impact  on  the  effi¬ 
ciency  gains  claimed  by  joined-wing  aircraft  [3].  He  estimates  that  joint  locations  from 
60%  to  100%  of  the  wing  span  should  be  evaluated  in  order  to  truly  find  an  optimum 
design.  Due  to  a  limitation  of  the  implementation  of  GEBT  used  in  this  research,  there 
is  no  consideration  of  a  joined  wing  with  any  structure  beyond  the  wing  joint  itself.  This 
research  considers  both  8  and  16  beam  element  models.  Regardless  of  the  number  of 
elements,  the  location  of  the  wing  tip  and  the  fore  and  aft  wing  root  locations  does  not 
change.  The  coordinates  for  the  wing  roots  are  shown  in  Table  3.1. 
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Tabic  3.1:  Simple  Joined- Wing  Model  Configuration 


Dimension 

Value 

(m) 

x(fore  root ) 

-7.173 

y(fore  root) 

-12.000 

z(fore  root ) 

-3.215 

x(aft  root ) 

7.173 

y(aft  root ) 

-12.000 

z(aft  root ) 

3.215 

c  chord 

1.500 

h  height 

0.150 

ts  skin  thickness 

0.015 

tw  web  thickness 

0.015 

The  model  under  examination  here  is  a  simple  box  beam  of  isotropic  material  as 
shown  in  Figure  3.2,  and  the  cross  sectional  elastic  constants  for  a  thin  walled  box  beam 
are 


1 

Tl 


T 

h 

1 


Figure  3.2:  Cross  section  of  simple  box  beam 
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EA  =  2  E[cts  +  ( h  —  2  ts)tw\ 


ch3 

(c  -  2 tw)(h  -  2 ts)3 

12 

12 

he3 

(h  -  2ts)(c  -  2 tw)3 

12 

12 

GJ  = 

4  GA2 

Jo  V* 

GJ  = 

2  Gc2h2 

( —  +  E) 

(3.1) 

(3.2) 

(3.3) 

(3.4) 

(3.5) 
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GA 


2  Get 


(3.6) 


XX 

GAZZ  =  2 G(h  -  2 ts)tw  (3.7) 

where  c  is  the  chord,  h  is  the  height,  ts  is  the  thickness  of  the  top  and  bottom  skins, 
and  tw  is  the  thickness  of  the  web.  For  the  simple  geometry  under  consideration,  Eqs. 
(3.1)  -  (3.7)  are  used  to  determine  the  cross-sectional  properties.  Equation  (3.4)  is  the 
general  form  of  the  torsion  constant  equation  for  thin  walled  sections  of  any  shape  [36]. 
Equation  (3.5)  is  the  form  used  in  the  calculations  in  this  research  as  applied  to  a  thin 
walled  box  beam.  For  more  complex  cross  sections,  VABS  is  well  suited  to  performing 
this  task  but  is  not  required  for  this  simple  cross  section. 

While  it  is  possible  to  apply  loads  to  any  node  in  the  model,  the  loads  are  limited 
here  to  static  loads  at  the  wing  tip  in  the  positive  global  Z-axis  direction.  A  small  500 
N  load  allows  comparison  of  displacement  results  from  linear  and  nonlinear  analyses  to 
ensure  correlation  early  in  the  setup  of  the  optimization  problem.  The  design  applied  tip 
load,  Pappiied  =  33,670  N,  is  that  described  by  Blair  in  Ref.  [10].  Table  3.2  shows  the  force 
and  moment  results  and  Table  3.3  shows  the  displacement  results  for  the  500  N  test  case 
for  both  the  nonlinear  solution  using  GEBT  and  a  linear  solution  using  Euler-Bernoulli 
beam  theory. 

The  element  forces  and  moments  are  in  the  local  beam  coordinate  system  while  the 
nodal  displacements  and  rotations  are  in  the  global  coordinate  system.  There  is  excellent 
agreement  between  these  two  sets  of  data.  This  research  focuses  on  the  displacement 
term  w3  at  node  5,  which  is  where  the  displacement  constraint  defined  in  Eq.(3.10)  is 
enforced.  The  nonlinear  solution  from  GEBT  is  0.001132  m.  When  the  equations  in 
GEBT  are  linearized,  the  tip  deflection  is  still  0.001132  m.  The  linear  solution  using 
Euler-Bernoulli  is  0.001141  m,  a  difference  of  0.8%.  The  linear  beam  model  used  is  not 
a  Timoshenko,  which  means  the  shear  correction  terms  are  not  considered.  The  cause  of 
this  difference  is  not  known  at  this  time,  although  the  nonlinear  and  linear  solutions  in 
GEBT  are  the  same,  so  linearization  of  GEBT  does  not  seem  to  be  the  cause. 

Larger  loads  are  expected  to  result  in  nonlinear  deformation  when  applied  to  the 
system.  Furthermore,  GEBT  should  yield  a  nonlinear  response  while  the  Euler-Bernoulli 
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Tabic  3.2:  500  N  test  load  force  and  moment  results  from  nonlinear  and  linear  analysis 

with  all  elements  set  to  initial  design  variable  values  as  given  in  Table  3.1 


Analysis 

Element 

Fi 

F2 

(N) 

Fs 

M1 

m2 

(N-m) 

m3 

1 

758 

-437 

-24.4 

128 

233 

-5450 

2 

758 

-437 

-24.4 

128 

145 

-3880 

3 

758 

-437 

-24.4 

128 

57 

-2310 

GEBT 

4 

758 

-437 

-24.4 

128 

-30 

-745 

5 

-758 

437 

-24.5 

-128 

30 

-745 

6 

-758 

437 

-24.4 

-128 

-57 

-2310 

7 

-758 

437 

-24.4 

-128 

-145 

-3880 

8 

-758 

437 

-24.4 

-128 

-232 

-5450 

1 

758 

-437 

-24.4 

128 

233 

-5450 

2 

759 

-437 

-24.4 

128 

145 

-3880 

3 

759 

-437 

-24.4 

128 

57 

-2310 

Euler- 

4 

759 

-437 

-24.4 

128 

-30 

-745 

Bernoulli 

5 

-759 

437 

-24.5 

-128 

30 

-745 

6 

-759 

437 

-24.4 

-128 

-57 

-2310 

7 

-759 

437 

-24.4 

-128 

-145 

-3880 

8 

-758 

437 

-24.4 

-128 

-232 

-5450 

solution  remains  linear  when  the  same  tip  load  is  applied  to  both  models.  Table  3.4  shows 
the  force  and  moment  results  while  Table  3.5  shows  the  displacement  and  rotation  results 
for  Papphed  —  33,670  N  for  both  the  nonlinear  solution  using  GEBT  and  a  linear  solution 
using  Euler-Bernoulli  beam  theory. 

Agreement  between  linear  and  nonlinear  analysis  at  Pappued  =  33,670  N  is  shown 
in  Table  3.5.  The  nonlinear  response  is  highlighted  best  by  examining  the  response  in 
the  U2  direction  and  rotation  62  at  node  5  in  Table  3.5.  Although  still  very  small,  the 
value  predicted  by  the  nonlinear  solver  for  d2  is  now  on  the  order  of  10-6,  an  increase  of 
10  orders  of  magnitude  over  the  numerically  zero  result  for  the  same  value  in  the  linear 
analysis.  Likewise,  the  112  deflection  has  appeared  and  is  the  same  order  of  magnitude  as 
the  other  nodes  in  the  nonlinear  case,  where  the  linear  solution  predicts  no  displacement 
in  the  global  Y  axis.  The  displacement  used  in  the  optimization,  w3  at  the  wing  tip  (node 
5),  is  76.26  mm  in  the  nonlinear  analysis  and  76.84  mm  in  the  linear  analysis,  a  0.8% 
difference.  Closer  inspection  of  the  data  shows  that  the  behavior  of  the  displacements  in 
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Tabic  3.3:  500  N  test  load  displacement  and  rotation  results  from  nonlinear  and  linear 

analysis  with  all  elements  set  to  initial  design  variable  values  as  given  in  Table  3.1 


Analysis 

Node 

Ml 

U2 

w 

u3 

0i 

e2 

(rad) 

o3 

2 

-42 

-2 

105 

5.86e-05 

-1.85e-05 

2.28e-05 

3 

-156 

-7 

381 

9.86e-05 

-2.46e-05 

4.09e-05 

4 

-317 

-8 

749 

1.20e-04 

-1.85e-05 

5.43e-05 

GEBT 

5 

-501 

0 

1132 

1.23e-04 

-4.99e-10 

6.31e-05 

6 

-317 

8 

749 

1.20e-04 

1.84e-05 

5.43e-05 

7 

-156 

7 

381 

9.87e-05 

2.46e-05 

4.09e-05 

8 

-42 

2 

105 

5.86e-05 

1.84e-05 

2.28e-05 

2 

-43 

-3 

110 

5.79e-05 

-1.83e-05 

2.29e-05 

3 

-158 

-8 

388 

9756e-05 

-2.44e-05 

4.10e-05 

4 

-320 

-9 

757 

1.19e-04 

-1.83e-05 

5.44e-05 

Euler- 

5 

-505 

0 

1141 

1.22e-04 

-1.70e-18 

6.31e-05 

Bernoulli 

6 

-320 

9 

757 

1.19e-04 

1.83e-05 

5.44e-05 

7 

-158 

8 

388 

9.75e-05 

2.44e-05 

4.10e-05 

8 

-43 

3 

110 

5.79e-05 

1.83e-05 

2.29e-05 

the  u2  and  u3  directions  as  well  as  F3  and  M\  show  a  noticeable  and  expected  departure 
from  linear  behavior  predicted  by  Eulcr-Bernoulli. 


3.2  Material  and  Geometric  Properties 

The  model  in  this  research  is  made  of  isotropic  prismatic  aluminum  beams.  The 
cross  section  is  a  thin-walled  box  with  four  variables  as  shown  in  Fig.  3.2.  Material 
properties  used  were  for  aluminum,  with  a  fixed  Young’s  modulus  E  of  70.0  GPa  and 
shear  modulus  G  of  26.31  GPa.  The  stress  criteria  used  were  von  Mises  stress,  defined 
as 


&VM 


(o 'll  -  cr22)2  +  (c22  -  CT 33)2  +  (cr33  -  CTn)2  +  6 (<jf2  +  g\3 


+  °3U 


(3.8) 


where  are  the  components  of  the  3D  stress  tensor.  Equation  (3.8)  can  be  applied 
to  both  the  nonlinear  and  linear  stresses  obtained  from  GEBT  and  Euler-Bernoulli, 
respectively.  For  the  stresses  computed  using  Eulcr-Bernoulli  linear  beam  theory,  the 
von  Mises  stresses  are  calculated  at  both  ends  of  each  beam  element  at  all  four  corners. 
These  points  are  shown  as  gray  triangles  in  Fig.  3.3.  In  the  Mathematica  implementation 
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Tabic  3.4:  33,670  N  load  force  and  moment  results  from  nonlinear  and  linear  analysis 

with  all  elements  set  to  initial  design  variable  values  as  given  in  Table  3.1 


Analysis 

Element 

Fi 

f2 

(kN) 

f3 

Mi 

m2 

TN-m) 

m3 

1 

51.1 

-29.5 

-1.81 

7.74 

16.0 

-366 

2 

51.1 

-29.4 

-1.61 

8.69 

9.46 

-261 

3 

51.2 

-29.4 

-1.49 

9.02 

3.64 

-155 

GEBT 

4 

51.2 

-29.3 

-1.44 

9.04 

-1.75 

-50.1 

5 

-50.9 

29.6 

-1.85 

-8.18 

2.30 

-50.2 

6 

-50.9 

29.6 

-1.80 

-8.21 

-4.12 

-156 

7 

-50.9 

29.5 

-1.67 

-8.56 

-10.1 

-262 

8 

-51.0 

29.4 

-1.48 

-9.52 

-1.53 

-368 

1 

51.1 

-29.4 

-1.62 

8.54 

15.5 

-367 

2 

51.1 

-29.5 

-1.63 

8.53 

9.65 

-262 

3 

51.1 

-29.5 

-1.63 

8.54 

3.81 

-156 

Euler- 

4 

51.1 

-29.6 

-1.63 

8.54 

-2.01 

-50.2 

Bernoulli 

5 

-51.1 

29.6 

-1.63 

-8.54 

2.01 

-50.2 

6 

-51.1 

29.5 

-1.63 

-8.54 

-3.81 

-156 

7 

-51.1 

29.5 

-1.63 

-8.54 

-9.65 

-262 

8 

-51.1 

29.4 

-1.62 

-8.54 

-15.5 

-367 

Table  3.5:  33,670  N  load  displacement  and  rotation  results  from  nonlinear  and  linear 

analysis  with  all  elements  set  to  initial  design  variable  values  as  given  in  Table  3.1 


Analysis 

Node 

Ml 

u2 

(mm) 

u3 

o2 

(rad) 

03 

2 

-2.85 

-0.18 

7.23 

3.94e-03 

-1.39e-03 

1.49e-03 

3 

-10.5 

-0.61 

25.9 

6.58e-03 

-1.75e-03 

2.72e-03 

4 

-21.3 

-0.80 

50.6 

8.02e-03 

-1.26e-03 

3.65e-03 

GEBT 

5 

-33.8 

-0.34 

76.3 

8.31e-03 

-2.27e-06 

4.25e-03 

6 

-21.3 

0.34 

50.3 

8.19e-03 

1.22e-03 

3.67e-03 

7 

-10.5 

0.32 

25.4 

6.71e-03 

1.56e-03 

2.79e-03 

8 

-2.86 

0.09 

6.97 

3.94e-03 

1.09e-03 

1.58e-03 

2 

-2.91 

-0.17 

7.39 

3.90e-03 

-1.23e-03 

1.54e-03 

3 

-10.6 

-0.51 

2.61 

6.57e-03 

-1.64e-03 

2.76e-03 

4 

-21.5 

-0.59 

51.0 

8.00e-03 

-1.23e-03 

3.67e-03 

Euler- 

5 

-34.0 

0.00 

76.8 

8.21e-03 

-1.21e-16 

4.25e-03 

Bernoulli 

6 

-21.5 

0.59 

51.0 

8.00e-03 

1.23e-03 

3.67e-03 

7 

-10.6 

0.51 

26.1 

6.57e-03 

1.64e-03 

2.76e-03 

8 

-2.91 

0.17 

7.39 

3.90e-03 

1.23e-03 

1.54e-03 
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of  GEBT  used  in  this  research,  stresses  are  not  returned  directly.  As  mentioned  earlier, 
the  GEBT  outputs  are  forces  and  moments.  These  forces  and  moments  are  calculated  at 
the  center  of  each  beam  element  in  the  beam  element’s  local  coordinate  system,  relieving 
the  user  of  the  need  to  apply  a  transformation  from  the  global  to  local  coordinate  system. 

Calculation  of  the  3D  nonlinear  stresses  and  strains  are  computed  by  the  current 
Fortran  implementation  of  VABS.  The  inputs  for  VABS  are  the  description  of  the  cross 
sectional  geometry,  material  properties,  and  forces  and  moments  applied  to  the  element. 
VABS  calculates  all  cross  sectional  properties  for  the  stiffness  matrix  prior  to  analysis, 
and  returns  3D  stress  and  strain  tensors  at  any  point  in  the  cross  section  desired  by 
the  user  at  the  midpoint  of  the  beam  element.  In  this  research,  a  cross  section  for  the 
thin  wall  box  beam  is  defined  using  20  nodes  which  are  grouped  into  four  elements. 
These  elements  comprise  the  top  and  bottom  skins  and  the  fore  and  aft  webs.  These 
elements  of  the  cross  section  are  shown  in  Fig.  3.4.  Each  element  is  defined  by  eight 
nodes,  six  of  which  are  shared  with  other  elements.  As  a  result,  the  stress  and  strain 
outputs  are  calculated  at  each  node  for  each  element,  resulting  in  32  outputs  for  stress 
recovery,  not  20.  The  stresses  reported  at  the  same  node  twice  are  averaged  to  give  a 
single  stress  tensor.  The  cross  section  is  shown  as  the  shaded  section  with  green  circles 
denoting  the  stress  recovery  points  used  in  the  von  Mises  stress  calculations  in  Fig. 
3.3.  The  corners  shown  at  the  midpoint  cross  section  relate  to  the  ends  of  the  beam 
for  comparison  purposes.  To  allow  comparison  between  the  GEBT/VABS  nonlinear 
and  the  Eulcr-Bernoulli  linear  stresses,  the  von  Mises  stresses  obtained  from  the  linear 
analysis  at  the  ends  of  the  beam  are  averaged.  This  approach  works  due  to  the  constant 
force  and  linear  moment  distributions  across  a  given  element.  These  averaged  linear  von 
Mises  stresses  compare  very  well  to  VABS  for  the  500  N  linear  test  case  at  the  midpoint 
cross-section  of  each  beam  element. 

3.3  Structural  Optimization  Problem 

For  each  beam  element  the  chord  and  height  of  the  wing  box  as  well  as  the  thickness 
of  the  skins  (upper  and  lower  surfaces)  and  webs  (fore  and  aft  surfaces)  are  all  design 
variables  as  shown  in  Fig.  3.2.  If  the  design  variables  of  the  eight  beams  listed  in  Fig.  1.2 
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Figure  3.3:  Simple  box  beam  showing  locations  of  von  Mises  stress  calculations 

are  linked,  there  are  only  four  design  variables.  In  the  case  where  each  of  the  eight  beam 
element’s  four  design  variables  are  independent  of  the  others,  there  are  thirty-two  design 
variables.  Upper  and  lower  bounds  were  defined  for  each  of  the  four  design  variables.  It 
is  important  to  note  that  the  constraint  on  the  height  of  the  box  beam  is  a  percentage  of 
the  current  value  of  the  chord.  As  a  result,  this  height  constraint  must  be  continuously 
updated  in  the  optimization  loop  as  values  of  the  chord  changes.  If  the  height  constraint 
is  not  updated,  it  is  possible  that  the  optimization  will  return  a  result  that  is  inconsistent 
with  this  constraint.  In  the  displacement  optimization  problem,  vertical  global  Z-axis 
deflection  at  the  wing  tip  is  the  design  constraint.  The  values  were  either  a  direct 
result  of  the  GEBT  nonlinear  analysis  or  Eq.  (3.11)  when  using  ESL.  Although  initially 
considered  as  two  separate  problems  during  development,  the  displacement  and  stress 
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Figure  3.4:  Cross  section  view  showing  nodes  used  by  VABS  to  calculate  stress  recovery 
points. 

constrained  optimization  problem  is  formulated  as: 


Find 
to  minimize 
subject  to 


b=  [c  h  ts  tw\T 

M(b ) 

|^Up|  ^  Hz— dir  allowed 
&  VM  <  O'  VM  allowed 

0.500  <  c  <  2.500  m 
0.05  c  <  h  <  0.15  c  m 
0.005  <  t,  <  0.015  m 


(3.9) 


0.005  <tw<  0.015  m 


The  mass  of  the  model  with  all  values  set  to  the  initial  design  variables  shown  in  Ta¬ 
ble  3.1  is  3764.8  kg.  A  MATLAB  sequential  quadratic  programming  (SQP)  optimization 
routine  is  used  for  the  optimization. 
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3-4  Integration  between  Mathematica  and  MATLAB 

Research  was  performed  in  a  Linux  environment  due  to  the  need  to  call  Mathe¬ 
matica  from  within  the  MATLAB  code  during  the  optimization  process.  A  script  used 
to  handle  writing  the  calculated  cross-sectional  properties  El,  GJ ,  et  al.  to  a  text 
hie.  Mathematica  is  then  called  using  a  shell  command  and  the  Mathematica  program 
(GEBT)  is  executed.  GEBT  initializes  its  variables,  then  reads  the  input  hies  created 
by  MATLAB  to  obtain  cross-sectional  data.  Nonlinear  analysis  is  then  conducted,  with 
displacement,  rotation,  force  and  moment  resultants  returned  in  a  single  matrix  that  is 
then  written  to  a  comma  separated  value  (CSV)  hie.  The  MATLAB  script  resumes  by 
reading  the  CSV  hie  and  transforming  the  Mathematica  output  into  MATLAB  format. 
The  values  returned  from  GEBT  are  at  the  mid  point  of  each  element,  so  a  linear  shape 
function  is  applied  to  determine  nodal  displacements  and  rotations  using  Eq.  (2.6) 


3.5  Nonlinear  Optimization  using  GEBT 

Two  different  approaches  to  the  optimization  problem  are  investigated.  In  the  hrst 
approach,  which  will  be  referred  to  as  GEBT,  a  full  nonlinear  response  optimization  is 
carried  out  using  GEBT  as  the  nonlinear  solver.  The  second  approach,  which  will  be 
referred  to  as  GEBT-ESL,  is  the  linear  response  optimization  using  ESL  in  conjunction 
with  the  GEBT  nonlinear  solver.  In  the  nonlinear  response  optimization,  the  optimiza¬ 
tion  is  driven  by  a  full  nonlinear  analysis.  The  optimization  starts  with  a  nonlinear 
analysis  of  the  model,  defines  upper  and  lower  bounds  on  the  design  variables,  and  then 
hands  the  model  off  to  SQP.  The  nonlinear  analysis  is  completed  by  a  Mathematica  im¬ 
plementation  of  GEBT  each  time  the  objective  and  constraint  equations  are  calculated 
within  SQP  using  the  current  values  of  all  design  variables.  The  joined-wing  analysis 
does  not  have  explicit  equations  programmed  that  can  be  used  to  calculate 


9  Ufip(li) 

9  bi 


(3.10) 


which  is  the  partial  derivative  of  the  wing  tip  displacement  with  respect  to  each  design 
variable.  The  gradient  based  search  method  employed  by  SQP  requires  a  sensitivity  anal- 
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ysis  with  respect  to  each  design  variable  be  performed  on  the  current  design,  including 
each  finite  difference  calculation  used  to  determine  the  search  direction. 

The  output  from  GEBT  is  then  processed  into  terms  that  are  used  in  the  constraint 
equations  for  displacement  and  stress  as  shown  in  Fig  3.5(a). 

3.6  Nonlinear  Optimization  using  GEBT  with  ESL 

The  second  approach  investigated  also  uses  the  Mathematica  implementation  of 
GEBT  for  nonlinear  analysis  in  an  outer  loop  to  determine  displacements  and  von  Mises 
stresses.  However,  a  set  of  equivalent  loads  calculated  with  a  linear  stiffness  matrix  before 
beginning  optimization  is  shown  in  Fig.  3.5(b)  and  is  used  in  a  linear  design  analysis. 
The  equivalent  loads  are  used  with  the  current  linear  stiffness  matrix  to  calculate  a  set 
of  approximate  nonlinear  displacements  from 

m  =  [jumru/y  (3.H) 

where  d  is  the  displacement  vector,  K^^b)  is  the  linear  stiffness  matrix  using  the  design 
variables  of  the  current  iteration,  b  is  the  vector  of  design  variables,  and  f*  is  the 
equivalent  static  loads  vector  for  displacements  that  was  calculated  in  the  outer  loop. 
From  the  displacements  d,  stresses  are  calculated  using  Euler-Bernoulli  beam  theory. 
The  element  stresses  are  then  used  to  calculate  the  von  Mises  stress  for  each  element 
using  Eq.  (3.8)  Within  the  optimization  inner  loop,  these  approximate  displacements 
and  stresses  will  be  used  in  the  evaluation  of  the  constraint  equations  in  Eq.  (3.10). 
Once  SQP  has  converged  to  a  solution,  the  results  are  returned  to  the  outer  loop,  where 
the  nonlinear  analysis  is  performed  again,  and  the  process  repeats  until  the  convergence 
criteria, 

£ converge 

is  met  in  the  outer  loop,  where  the  ||/g  ||^  term  is  the  infinity  norm  of  the  equivalent 
static  load  vector  for  the  current  nonlinear  evaluation  and  H/eJI^1  is  the  infinity  norm 
of  equivalent  static  load  vector  for  the  previous  nonlinear  evaluation.  This  convergence 
tolerance  is  satisfied  as  shown  in  Fig.  3.5. 
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(a)  GEBT 


(b)  GEBT-ESL 


Figure  3.5:  Flow  charts  showing  the  difference  between  pure  nonlinear  optimization 

(GEBT)  and  nonlinear  optimization  with  equivalent  static  loads  (GEBT-ESL). 
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3.7  Limit  Load  Calculations 


Initial  efforts  to  perform  the  nonlinear  optimization  revealed  peculiarities  in  the 
joined-wing  optimization  problem.  In  the  first  run  attempted,  the  optimizer  set  all  design 
variables  to  minimum  gage  after  the  first  step  in  an  attempt  to  minimize  the  weight  of 
the  structure.  The  problem  arises  that,  when  the  minimum  gage  design  as  defined  in  this 
research  is  sent  to  GEBT  for  a  nonlinear  response  analysis,  the  analysis  fails.  The  issue 
is  hinted  at  in  a  message  returned  by  Mathematica,  stating  that  the  Jacobian  used  in 
the  calculation  of  the  solution  is  singular,  and  that  the  design  should  be  adjusted  in  an 
effort  to  obtain  a  well  behaved  Jacobian.  After  failing  to  harness  a  flag  from  this  error 
message  as  a  means  of  triggering  a  ‘bad  design’  indicator  in  MATLAB,  an  alternative 
approach  is  presented. 


Figure  3.6:  Load-deflection  curve  showing  linear  and  nonlinear  behavior.  PcriUcai  is 

located  at  the  peak  of  the  nonlinear  curve  and  indicates  failure  of  the  design. 

Elasticity  theory  describes  the  load-deformation  curve  and  uses  it  to  show  linear- 
elastic  and  nonlinear-elastic  regions  of  deformation.  If  the  load  on  a  specimen  continues 
to  increase  beyond  the  yield  stress,  behavior  of  the  material  cannot  be  predicted  with 
GEBT.  The  analysis  by  GEBT  shows  this  failure  in  the  form  of  the  poorly  conditioned 


stiffness  matrix.  The  stiffness  matrix  becomes  poorly  conditioned  due  to  the  geometric 
nonlinearity  effects.  Geometric  nonlinearities  con  manifest  themselves  either  in  an  un¬ 
stable  stiffness  matrix  due  to  buckling  (cannot  be  inverted)  or  by  the  presence  of  large 
nonlinear  strain  terms  which  render  the  stiffness  matrix  numerically  ill  conditioned  (re¬ 
sults  unreliable).  The  solution  to  determining  this  failure  criterion  is  to  generate  the 
load-deformation  curve  each  time  the  constraints  are  evaluated  to  find  the  value  of  a 
PcHticai  that  is  defined  as  the  point  where  the  deflection  returned  from  the  nonlinear 
solver  is  less  than  the  previous  value  on  the  curve.  A  graphical  example  of  this  location 
is  provided  in  Fig.3.6  where  the  nonlinear  curve  is  at  a  maximum  value  of  approximately 
6.3  m.  Furthermore,  it  is  clear  that  operating  near  this  region  of  the  curve  could  be 
disastrous  in  a  real  world  design,  so  a  safety  factor  Psafe  is  defined 

P 

Pi  criticcil  / o  i  o\ 

safe  =  15  {0.16) 

and  used  as  a  constraint  on  the  optimization  problem  as 

Papplied  Psafe •  (3.14) 

It  is  possible  to  calculate  a  nonlinear  optimization  of  this  problem  with  less  com¬ 
putational  effort.  Using  rounded  off  values  of  the  design  variables,  a  comparison  can  be 
made  that  eliminates  a  large  percentage  of  the  nonlinear  analysis  overhead.  For  example, 
the  following  set  of  design  variables  are  handed  to  the  optimizer  for  step  kn : 

c  =  1.2573874 

h  =  0.1088356  (3.15) 

ts  =  0.0623400 
tw  =  0.0500000 
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and  are  rounded  off  to  four  decimal  places.  The  rounded  off  values  for  step  kn  are 

c  =  1.2574 

h  =  0.1088  (3.16) 

ts  =  0.0623 
tw  =  0.0500 

which  are  stored  for  comparison  during  the  next  step  in  the  optimization.  During  a 
finite  difference  method,  the  small  value  added  to  the  design  variables  is  on  the  order  of 
10~6  or  smaller.  In  this  example,  10-6  is  used,  so  that  the  design  variables  for  the  next 
iteration  kn+\  are 

c  =  1.2573884 

h  =  0.1088366  (3.17) 

t9  =  0.0623410 
tw  =  0.0500010 

which  rounded  off  to  four  decimal  places  are  the  same  values  shown  in  Eqs.  (3.17). 
When  a  comparison  is  made  between  the  rounded  set  of  design  variables  for  step  kn 
and  kn+\  they  are  found  to  be  equivalent  regarding  the  change  in  the  shape  of  the  load- 
displacement  curve.  There  is  still  sufficient  change  when  using  the  unrounded  design 
variables  to  calculate  the  required  displacement  and  stress  sensitivities  needed  for  SQP 
to  choose  the  appropriate  search  direction.  Now  that  the  load-deflection  curve  is  not 
changing  with  the  finite  difference  calculations,  it  is  then  extended  to  the  calculation  of 
the  critical  and  safe  loads  for  a  given  design.  Since  the  results  for  Pcriticai  and  Psafe  do 
not  change  after  the  first  finite  difference  calculation,  these  nonlinear  analyses  do  not 
need  to  be  performed  and  can  be  removed  while  using  the  previous  values  for  Pcriticai 
and  Psafe  in  a  constraint. 

Upon  further  investigation,  it  is  clear  that  knowing  the  actual  limit  load  is  unnec¬ 
essary  if  the  design  is  safe.  This  is  because  the  direction  of  search  taken  in  the  optimizer 
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will  not  be  based  on  a  constraint  violation  of  the  limit  load.  It  is  enough  to  know  that 
Psafe  is  not  less  than  the  applied  load.  With  this  in  mind,  further  time  savings  can  be 
had  by  limiting  the  upper  limit  of  the  load-deformation  curve  calculation  to  no  more  than 
1.5  times  the  applied  tip  load  and  implementing  an  active  set  strategy  for  the  gradient 
calculation.  This  will  ensure  adequate  determination  of  the  constraint  violations  while 
removing  all  extraneous  nonlinear  calls  in  a  nonlinear  response  optimization. 

3.8  Limit  Load  Estimation  Techniques 

3.8.1  Polynomial  Curve  Fit  with  Equivalent  Static  Loads  at  Limit  Load.  In 
order  to  make  a  fair  assessment  in  the  gains  offered  by  the  use  of  ESL  a  means  of  including 
the  limit  load  in  the  GEBT-ESL  inner  loop  is  desired.  Also,  inclusion  of  the  limit  load 
constraint  should  also  guarantee  that  the  limit  load  will  not  be  violated  by  the  final 
design  returned  from  the  optimizer.  A  second  order  polynomial  curve  fit  is  an  important 
piece  to  an  inner  loop  approximation  of  the  limit  load  that  can  be  applied  as  a  constraint. 
Instead  of  defining  the  limit  load  curve  with  a  constant  step  size,  a  selection  of  four  to 
five  points  in  the  vicinity  of  the  applied  load  and  desired  limit  load  are  evaluated  using 
GEBT.  These  loads  then  return  a  displacement  vector  for  the  loads  under  consideration. 
The  use  of  a  quadratic  curve  fit  algorithm  will  result  in  the  polynomial 

P(u)  =  a0  +  aiu  +  a-2u2  (3.18) 

where  a o,  cq,  and  cq  are  the  coefficients  when  the  curve  fit  algorithm  is  supplied  with  the 
calculated  displacements.  This  means  Eq.  (3.18)  will  return  a  load  when  supplied  with  a 
new  displacement.  To  determine  the  displacement  in  the  inner  loop,  a  separate  set  of  ESL 
are  calculated  in  the  outer  loop  for  the  limit  load  calculation  using  the  maximum  load 
used  in  the  calculation  of  the  curve  fit  data  points  in  GEBT  to  return  the  true  nonlinear 
wing  tip  displacement.  Differentiating  Eq.  (3.18)  with  respect  to  tip  displacement, 

dP(u )  ,  . 

— - —  =  a i  +  2  *  cqw  (3.19) 


31 


which  if  set  equal  to  zero  can  be  solved  for  Uiirnit.  Then  using  this  new  value  for  Uiirnit 
in  Eq.  (3.18)  will  produce  an  estimate  of  Pumit  as  shown  in  Fig.  3.7.  Using  uumit 
an  ESL  can  be  calculated  for  the  limit  load  displacement.  Within  the  inner  loop,  the 
ESL  at  the  limit  load  will  generate  different  deflections  at  the  limit  load  as  the  design 
variables  are  updated.  Results  with  this  Purmt  estimation  technique  within  the  GEBT- 
ESL  displacement  constraint  problem  are  shown  in  Table  4.4 

The  use  of  this  method  of  Pumit  estimation  does  yield  some  improvements  over  the 
calculation  where  Pumit  is  ignored  as  a  constraint.  In  stiffer  structures,  the  constraint 
appears  to  allow  for  a  more  efficient  design,  but  at  the  cost  of  additional  computational 
cost  as  shown  in  Table  4.4.  This  tradeoff  of  computational  effort  for  much  improved 
reliability  is  considered  worth  the  investment.  By  comparing  Fig.  4.1  to  Fig.  4.2  it  is 
clear  there  is  improvement  in  the  optimized  mass  of  the  structure  when  the  wing  tip 
deflection  is  constrained  to  0.25  and  0.35  m.  At  the  other  end  of  the  scale,  use  of  the 
limit  load  estimation  led  to  results  at  displacement  constraint  values  that  were  previously 
unattainable  for  the  32-DV  and  64-DV  problems. 


3.8.2  Polynomial  Curve  Fit  with  Updated  Slope  Term.  Another  possible  ap¬ 
proach  to  the  estimation  of  Pumit  is  to  update  the  Gq  or  slope  term  in  Eq.  (3.18)  inside 
the  inner  loop.  The  setup  is  the  same  in  the  outer  loop  as  described  in  the  previous 
section.  When  the  polynomial  is  handed  to  the  inner  loop,  a  different  calculation  is 
performed.  A  ratio 


~  ^ applied 

d\  =  — — - d\ 

^ applied 


(3.20) 


where  uappiied  is  the  tip  deflection  calculated  by  GEBT  in  the  outer  loop  at  the  applied 
load,  cii  is  the  coefficient  of  the  polynomial  calculated  in  the  outer  loop,  uappued  is  the 
inner  loop  wing  tip  deflection  estimated  using  ESL  at  the  applied  load,  and  cq  is  the 
updated  value  for  the  slope  in  the  inner  loop.  The  change  in  slope  allows  an  update  to 
the  limit  load  in  the  linear  analysis.  This  update  method  yields  results  similar  to  those 
provided  in  the  previous  section  with  no  obvious  advantage  at  this  time. 
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Figure  3.7:  Plot  showing  the  nonlinear  load  deflection  curve  using  the  incremental 

step  method  (blue  solid  line),  the  linear  load  deflection  curve  (green  solid  line),  and  the 
2nd  order  polynomial  curve  fit  for  the  limit  load  using  a  sampling  of  5  points. 


3.9  Trust  Region  Strategy 

Move  limits  play  a  vital  role  in  the  ability  of  the  nonlinear  response  optimization 
to  converge.  The  fixed  reduction  of  the  move  limit  does  not  always  allow  the  overall 
MATLAB  routine  to  converge  to  a  solution  that  improved  upon  the  4-DV  case  when 
starting  with  the  same  initial  design  variables,  even  though  SQP  successfully  converged  in 
the  inner  optimization  loop.  Dynamic  move  limits  can  be  used  to  allow  the  optimization 
routine  to  adjust  the  step  size  while  searching  for  the  optimal  solution,  but  require 
additional  information  to  work  properly.  A  definition  of  when  to  expand,  contract,  or 
keep  the  current  move  limits  is  needed.  Introduction  of  a  trust  region  strategy  proved  to 
be  an  acceptable  solution  where  the  trust  ratio  (ft  is  defined  as 

iftk- 1  -  'iftk 
(ft  =  - — 

iftk- i  -  iftk 


(3.21) 


33 


where  V’fc-i  is  the  value  of  merit  function  based  on  the  nonlinear  analysis  before  the 
optimization  loop,  0*,  is  the  value  of  the  merit  function  based  on  the  nonlinear  analysis 
after  the  optimization  loop,  and  0*,  is  the  value  of  the  merit  function  using  the  equivalent 
static  loads.  The  merit  function  used  in  the  displacement  constraint  problem  was  defined 
as 

0  \J U tip  allow )  (3.22) 

where  utip  is  the  vertical  displacement  of  the  wing  tip  node  and  utip  allow  is  the  allowed 
displacement  of  the  same  node.  The  trust  region  strategy  allows  the  outer  loop  of 
the  optimization  to  dynamically  choose  whether  to  maintain,  expand,  or  contract  the 
allowable  move  limits  for  the  inner  loop  ESL  optimization  based  on  the  value  of  0  after 
each  result  from  the  inner  loop  is  returned.  If  a  set  of  design  variables  are  rejected,  the 
values  used  to  enter  the  inner  optimization  loop  are  restored  and  optimized  again  with 
contracted  move  limits.  If  0  indicates  an  improvement  in  the  design,  the  move  limits  are 
kept  the  same.  If  0  shows  a  great  improvement  in  the  design,  then  the  move  limits  are 
expanded  in  an  attempt  to  allow  the  design  to  reach  the  optimum  more  quickly. 
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IV.  Results  and  Discussion 


THE  results  for  the  nonlinear  optimization  problem  are  presented  here.  Wing  box 
shape  and  sizing  is  accomplished  using  the  methods  described  in  Chapter  3  and 
integration  of  GEBT  and  ESL  is  an  effective  combination  to  reduce  computational  ef¬ 
fort  while  maintaining  accuracy  in  the  final  solution.  There  are  also  some  benefits  to 
increasing  the  number  of  beam  elements  in  the  structure,  at  the  expense  of  increasing 
the  computational  expense  compared  to  smaller  numbers  of  beam  elements/degrees  of 
freedom.  Finally,  although  there  is  no  mention  in  the  literature  regarding  limit  load 
criteria  or  constraints,  it  is  shown  here  to  be  an  important  safety  consideration  which 
must  be  accounted  for  in  the  design  optimization  process. 

4-1  Four  Design  Variables  -  Displacement  Constraints 

4-1.1  Nonlinear  Optimization  (GEBT).  The  full  nonlinear  response  optimiza¬ 
tion  solution  with  uzauow  =  0.25  m  is  shown  in  Table  4.1.  The  mass  of  the  optimized 
structure  is  1183.1  kg.  Note  that  the  results  for  all  elements  are  the  same  due  to  the 
design  variable  linking  in  the  4-DV  case.  The  results  presented  in  Table  4.1  are  the 
same  when  the  limit  load  is  calculated  every  time  the  constraints  are  evaluated  and 
takes  4  hours  24  minutes  to  complete.  Eliminating  unneccessary  calculation  of  the  full 
load-deflection  curve  from  the  finite  differencing  calls  using  the  rounded  design  variable 
comparison  resulted  in  a  run  time  of  33  minutes  34  seconds,  a  92%  time  savings  over  the 
full  analysis.  A  comparison  of  different  approaches  to  the  4-DV  dispacment  constraint 
optimization  problem  is  shown  in  Table  4.2. 

Computationally,  the  addition  of  the  limit  load  constraint  is  significant.  The  limit 
load  for  the  initial  set  of  design  variables  is  more  than  250,000  N,  which  requires  50  calls 
to  the  nonlinear  solver  with  a  step  size  of  5,000  N  to  determine  this  value.  Other  step 
sizes  will  change  the  required  number  of  calls,  but  as  Psafe  approaches  the  applied  load,  a 
finer  step  size  may  be  required  to  accurately  determine  if  the  design  meets  the  limit  load 
constraint.  Computing  the  limit  load  calculation  for  every  step  in  the  finite  differencing 
routine  dramatically  increases  the  time  required  to  perform  the  optimization,  on  the 
order  of  hours  instead  of  minutes  or  seconds. 


35 


Table  4.1:  Solution  to  4-DV  displacement  constraint  nonlinear  response  optimization 
where  uzaUow  =  0.25  m 


Element 

c 

(m) 

h 

(m) 

^ skin 

(m) 

t"web 

(m) 

1 

1.3339 

0.2005 

0.0050 

0.0050 

2 

1.3339 

0.2005 

0.0050 

0.0050 

3 

1.3339 

0.2005 

0.0050 

0.0050 

4 

1.3339 

0.2005 

0.0050 

0.0050 

5 

1.3339 

0.2005 

0.0050 

0.0050 

6 

1.3339 

0.2005 

0.0050 

0.0050 

7 

1.3339 

0.2005 

0.0050 

0.0050 

8 

1.3339 

0.2005 

0.0050 

0.0050 

4-1-2  Nonlinear  Optimization  with  Linear  Approximation  (GEBT-ESL).  When 
equivalent  static  loads  are  introduced,  they  replace  the  full  nonlinear  analysis  inside  SQP 
with  a  linear  approximation.  This  means  inside  the  inner  loop,  there  is  no  difficulty  with 
a  cumbersome  constraint  evaluation.  As  a  result  of  using  ESL  and  not  requiring  a  limit 
load  constraint,  the  time  savings  is  three  orders  of  magnitude  between  the  full  nonlinear 
GEBT  and  the  GEBT-ESL  solutions.  It  should  be  noted  that  the  optimization  for 
GBT-ESL  did  not  evaluate  the  limit  load  constraint  in  any  way.  This  frees  GEBT-ESL 
from  the  computational  effort  required  for  nonlinear  optimization  that  is  largely  due  to 
calculating  this  constraint.  Also,  it  is  shown  in  Table  4.1  and  Table  4.2  that  the  final 
solutions  for  each  method  return  the  same  objective  function  value,  Psafe,  and  design 
variable  values.  The  only  changes  are  in  computational  effort  and  time  required.  Note 
that  the  values  for  Psafe  shown  are  calculated  after  the  optimization  is  complete  using 
the  final  design  variable  vector  b  and  is  not  included  in  the  time  required  to  arrive  at  the 
solution.  The  results  for  the  4-DV  displacement  constraint  problem  using  GEBT-ESL 
without  limit  load  calculations  are  shown  in  Table  4.3  and  Fig.  4.1. 

4-2  Thirty-two  Design  Variables  -  Displacement  Constraints 

4-2.1  Move  Limits  and  Move  Reduction.  Increasing  the  number  of  design 
variables  for  the  displacement  constraint  scenario  is  studied.  In  this  case,  eight  elements 
with  four  design  variables  each  are  all  independent  of  each  other  for  a  total  of  thirty-two 
design  variables  (32-DV).  Initially,  it  was  believed  that  the  code  would  be  identical  to  that 
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Figure  4.1:  Pareto  curve  showing  change  in  optimized  mass  as  the  tip  deflection  con¬ 
straint  is  varied  for  4,  32,  and  64-DV  cases 


Figure  4.2:  Pareto  curve  showing  change  in  optimized  mass  as  the  tip  deflection  con¬ 
straint  is  varied  for  4,  32,  and  64-DV  cases  with  a  2nd  order  polynomial  curve  fit  estimate 
of  Pumit  applied  as  a  constraint  in  the  inner  loop. 
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Tabic  4.2:  Comparison  of  different  solution  techniques  when  applied  to  the  four  design 
variable  optimization  problem  with  a  displacement  constraint. 


Solution 

Method 

^ zallow 

(m) 

fobj 

(kg) 

Move 

Limit 

Move 

Reduction 

Time 

Required 

Nonlinear 

Calls 

Psafe 

(kN) 

GEBT 

0.25 

1183.1 

0.09 

0.7 

4  h  23  m 

>10,000 

297 

GEBT-DV  compare 

0.25 

1183.1 

0.09 

0.7 

33  m  34  s 

1821 

297 

GEBT-DV  compare 
with  2  x  Load  check 

0.25 

1183.1 

0.09 

0.7 

8  40  s 

304 

297 

GEBT-ESL 

0.25 

1183.1 

0.09 

0.7 

8  s 

4 

297 

of  the  4-DV  problem  with  independent  values  for  c,  h,  ts,  and  tw.  When  an  optimization 
is  performed  without  move  limits  or  use  of  the  limit  load  calculation  described  in  Section 
3.7,  the  design  shown  in  Table  4.5  is  returned,  with  fabj  =  861.6  kg  and  Psafe  —  20  kN. 
Time  required  for  the  analysis  is  1  minutes  43  seconds,  which  is  reasonable  considering 
the  increased  number  of  design  variables  over  the  4-DV  case.  The  critical  issue  here  is 
the  violation  of  Psafe  when  checked  at  the  end  of  the  optimization.  The  advantage  in 
computational  time  seen  in  the  4-DV  case  in  GEBT-ESL  disappears  with  the  increase 
in  the  number  of  design  variables. 

Different  approaches  could  potentially  solve  this  issue,  however  caution  must  be 
exercised  in  selecting  the  method  lest  the  advantages  gained  in  using  the  equivalent 
static  loads  will  disappear.  Most  obviously,  the  use  of  the  limit  load  calculation  originally 
described  in  the  4-DV  nonlinear  analysis  could  be  used.  It  is  still  very  costly,  however,  to 
calculate  the  load-deformation  curve  for  each  step  of  the  optimization.  This  approach  is 
not  recommended  as  a  result.  A  second  possibility  is  the  use  of  a  trust  region  strategy  in 
the  outer  loop  of  the  ESL  routine  that  will  detect  a  violation  of  the  limit  load  constraint 
and  reduce  the  allowable  change  in  design  variables  for  a  given  nonlinear  analysis.  The 
reduction  of  the  move  limit  continues  until  either  the  next  step  does  not  violate  the  limit 
load  or  the  solution  has  converged.  A  third  option  is  that  move  limits  could  be  prescribed 
for  each  individual  design  variable,  or  made  relative  to  each  design  variable  as  opposed 
to  using  absolute  move  limits  for  all  DVs.  If  using  relative  move  limits  and  the  initial 
value  of  c  is  1.5  m,  while  ts  is  0.015  m,  it  may  be  desirable  to  allow  the  value  of  c  to 
vary  ±0.1  m  while  ts  is  only  allowed  to  vary  ±  0.001  m  per  iteration.  This  technique 
becomes  more  important  the  larger  the  difference  in  magnitude  of  the  design  variables 
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Table  4.3:  Optimization  results  for  beam  elements  with  displacement  constraints. 


Solution 

—  allow 

fobj 

Move 

Move 

Time 

Nonlinear 

Psafe 

Method 

(m) 

(kg) 

Limit 

Reduction 

Required 

Calls 

(kN) 

0.25 

1183.1 

0.1 

0.5 

12  s 

4 

297 

0.35 

1057.1 

0.1 

0.5 

9  s 

4 

207 

0.45 

972.3 

0.1 

0.5 

10  s 

5 

160 

0.55 

909.8 

0.1 

0.5 

12  s 

6 

127 

4-DV 

0.65 

861.1 

0.1 

0.5 

12  s 

6 

107 

GEBT-ESL 

0.75 

821.7 

0.1 

0.5 

13  s 

7 

93 

0.85 

789.0 

0.1 

0.5 

14  s 

7 

80 

0.95 

761.3 

0.1 

0.5 

14  s 

7 

73 

1.05 

737.5 

0.1 

0.5 

14  s 

7 

70 

1.15 

716.7 

0.1 

0.5 

17  s 

7 

63 

1.25 

698.4 

0.1 

0.5 

19  s 

9 

57 

0.25 

1177.8 

0.09 

0.7 

55  s 

11 

120 

0.35 

804.9 

0.09 

0.7 

77  s 

14 

127 

0.45 

727.7 

0.10 

0.5 

30  s 

7 

97 

0.55 

691.6 

0.09 

0.7 

58  s 

11 

80 

32-DV 

0.65 

664.4 

0.10 

0.5 

73  s 

18 

77 

GEBT-ESL 

0.75 

642.0 

0.08 

0.5 

30  s 

7 

60 

0.85 

623.7 

0.09 

0.7 

40  s 

8 

53 

0.95 

608.6 

0.10 

0.5 

25  s 

7 

47 

1.05 

595.2 

0.09 

0.7 

38  s 

8 

43 

1.15 

584.3 

0.10 

0.5 

58  s 

11 

40 

1.25 

- 

- 

- 

- 

- 

- 

0.25 

843.6 

0.07 

0.5 

210  s 

6 

147 

0.35 

771.6 

0.07 

0.5 

162  s 

6 

107 

0.45 

723.5 

0.07 

0.5 

124  s 

6 

87 

0.55 

700.4 

0.07 

0.5 

254  s 

8 

60 

64-DV 

0.65 

663.1 

0.07 

0.5 

164  s 

7 

53 

GEBT-ESL 

0.75 

641.9 

0.07 

0.7 

161  s 

7 

47 

0.85 

625.0 

0.08 

0.5 

265  s 

10 

40 

0.95 

610.4 

0.08 

0.5 

143  s 

7 

40 

1.05 

598.5 

0.09 

0.7 

141  s 

7 

36 

39 


Tabic  4.4:  Optimization  results  for  beam  elements  with  displacement  constraints  and 
Psafe  polynomial  curve  fit  estimation  routine  incorporated 


Solution 

—  allow 

fobj 

Move 

Move 

Time 

Nonlinear 

Psafe 

Method 

(m) 

(kg) 

Limit 

Reduction 

Required 

Calls 

(kN) 

0.25 

1183.1 

0.1 

0.5 

21  s 

16 

293 

0.35 

1057.1 

0.1 

0.5 

21  s 

16 

207 

0.45 

972.3 

0.1 

0.5 

21  s 

16 

157 

0.55 

909.8 

0.1 

0.5 

21  s 

16 

133 

4-DV 

0.65 

861.1 

0.1 

0.5 

21  s 

16 

110 

GEBT-ESL 

0.75 

821.7 

0.1 

0.5 

27  s 

21 

97 

0.85 

789.0 

0.1 

0.5 

27  s 

21 

87 

0.95 

761.4 

0.1 

0.5 

28  s 

21 

73 

1.05 

737.5 

0.1 

0.5 

27  s 

21 

67 

1.15 

716.9 

0.1 

0.5 

27  s 

21 

60 

1.25 

698.4 

0.1 

0.5 

27  s 

21 

53 

0.25 

985.1 

0.09 

0.7 

310  s 

178 

120 

0.35 

792.6 

0.09 

0.7 

194  s 

126 

127 

0.45 

727.6 

0.10 

0.5 

68  s 

42 

97 

0.55 

691.7 

0.09 

0.7 

58  s 

35 

80 

32-DV 

0.65 

663.8 

0.10 

0.5 

198  s 

133 

77 

GEBT-ESL 

0.75 

640.9 

0.08 

0.5 

207  s 

133 

60 

0.85 

622.0 

0.09 

0.7 

201  s 

133 

53 

0.95 

606.0 

0.10 

0.5 

243  s 

133 

47 

1.05 

596.7 

0.09 

0.7 

228  s 

140 

43 

1.15 

584.0 

0.10 

0.5 

198  s 

109 

40 

1.25 

576.0 

0.10 

0.5 

348  s 

231 

40 

0.25 

864.0 

0.9 

0.7 

869  s 

91 

47 

0.35 

772.3 

0.4 

0.7 

162  s 

28 

95 

0.45 

723.8 

0.4 

0.7 

183  s 

28 

69 

0.55 

689.9 

0.4 

0.7 

196  s 

35 

69 

64-DV 

0.65 

663.0 

0.5 

0.5 

223  s 

42 

59 

GEBT-ESL 

0.75 

642.0 

0.5 

0.7 

350  s 

63 

48 

0.85 

624.9 

0.5 

0.5 

172  s 

35 

45 

0.95 

610.9 

0.5 

0.5 

171  s 

35 

41 

1.05 

598.0 

0.5 

0.5 

231  s 

49 

37 

1.15 

598.0 

0.5 

0.8 

318  s 

49 

37 

40 


being  considered,  as  it  will  allow  slower  changes  to  the  design  that  will  more  likely  keep 
the  design  in  the  feasible  design  space  as  iterations  occur. 

Using  move  limits  while  ignoring  inner  loop  calculation  or  approximation  of  the 
limit  load,  GEBT-ESL  returns  designs  that  meet  the  Pcriucai  value  but  still  violate  the 
Psafe  constraint. There  are  also  other  difficulties  that  arise.  When  using  the  initial  design 
variable  values  in  Table  3.1,  and  an  outer  loop  move  limit  on  all  design  variables  of  0.15 
m,  the  design  shown  in  Table  4.6  is  the  result.  Of  greatest  concern  is  the  fact  that 
neither  the  front  or  aft  wing  tapers  smoothly  from  the  root  to  the  tip  in  the  chord  and 
height  variables.  Also,  the  objective  function  value  fQbj  =  1706.3  kg  is  approximately 
a  500  kg  increase  over  the  4-DV  solution  already  calculated,  so  this  cannot  be  the  true 
global  minimum.  The  fixed  move  limits  do  not  afford  the  optimization  problem  the  time 
needed  to  converge  to  a  true  minimum  solution.  When  the  move  limits  were  expanded 
to  allow  SQP  to  move  the  design  more  quickly,  the  final  result  did  not  meet  the  limit 
load  constraint. 

Experimentation  with  different  starting  design  variables  showed  that  it  is  possible 
to  improve  upon  the  4-DV  GEBT-ESL  results.  When  the  initial  values  for  the  opti¬ 
mization  were  set  to  the  values  obtained  for  the  4-DV  solution,  GEBT-ESL  was  able  to 
improve  upon  the  4-DV  result.  For  move  limit  =  0.1  m  and  move  reduction  =  70  %,  fQbj 
=  1097.5  kg  and  had  a  Psafe  —  37  kN.  This  is  a  valid  design  satisfying  all  constraints, 
even  though  no  limit  load  calculations  were  being  performed  during  the  optimization. 
As  a  result,  only  1  min  43  sec  were  required  to  converge  to  the  solution. 

Note  the  only  issue  in  Table  4.7  is  that  the  the  value  of  h  for  element  1  (0.1340  m)  is 
less  than  the  value  for  element  2  (0.1780  m).  This  violates  the  smooth  taper  from  root  to 
tip  sought  after  in  the  design.  While  not  strictly  necessary,  intuition  would  indicate  that 
a  smooth  taper  would  be  the  lightest  design  in  both  a  displacement  or  stress  constrained 
optimization  problem. 

Adjusting  the  move  limit  to  0.05  m  and  leaving  the  move  reduction  at  70%  returns 
the  best  solution  seen  to  this  point.  The  design  satisfies  all  constraints,  has  a  smooth 
taper  from  root  to  tip  in  both  wings,  and  does  not  require  calculation  of  the  limit  load 
within  the  algorithm.  The  value  of  the  objective  function  is  fQbj  =  1122.9  kg,  Psafe  —  40 


41 


Tabic  4.5:  Design  variables  by  element  for  32-DV  GEBT-ESL  design  with  no  move 
limits  or  limit  load  calculation. 


Element 

chord 

height 

skin 

web 

1 

1.8943 

0.2842 

0.0050 

0.0083 

2 

1.6624 

0.2494 

0.0050 

0.0082 

3 

1.3128 

0.1969 

0.0050 

0.0082 

4 

0.7047 

0.1057 

0.0050 

0.0084 

5 

0.5000 

0.0250 

0.0050 

0.0010 

6 

0.5000 

0.0250 

0.0050 

0.0010 

7 

0.5000 

0.0250 

0.0050 

0.0010 

8 

0.5000 

0.0250 

0.0050 

0.0010 

Table  4.6:  Design  variables  by  element  for  32-DV  GEBT-ESL  design  with  move  limit 
=  0.15  m  and  move  reduction  =  70  %  with  no  limit  load  calculation.  Initial  design 
variables  are  those  used  in  Table  3.1. 


Element 

chord 

height 

skin 

web 

1 

1.3982 

0.1332 

0.0089 

0.0177 

2 

1.4477 

0.1487 

0.0089 

0.0145 

3 

1.0647 

0.1439 

0.0089 

0.0136 

4 

0.8934 

0.0899 

0.0089 

0.0090 

5 

0.8883 

0.0888 

0.0089 

0.0089 

6 

0.9214 

0.1136 

0.0089 

0.0113 

7 

1.0496 

0.1493 

0.0089 

0.0114 

8 

0.9901 

0.1485 

0.0089 

0.0121 
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Tabic  4.7:  Design  variables  by  element  for  GEBT-ESL  design  with  move  limit  =  0.10 
m  and  move  reduction  =  70  %  with  no  limit  load  calculation.  Initial  design  variables  are 
the  solution  to  the  4-DV  0.25  m  tip  displacement  constraint  optimization  problem. 


Element 

chord 

height 

skin 

web 

1 

1.7252 

0.1340 

0.0050 

0.0113 

2 

1.1904 

0.1780 

0.0050 

0.0077 

3 

0.8935 

0.1340 

0.0050 

0.0062 

4 

0.8935 

0.1340 

0.0050 

0.0058 

5 

0.8935 

0.1340 

0.0050 

0.0058 

6 

1.0401 

0.1534 

0.0050 

0.0073 

7 

1.3493 

0.1942 

0.0050 

0.0088 

8 

1.4743 

0.2151 

0.0050 

0.0087 

kN,  and  time  required  is  1  min  22  sec.  Initial  design  variables  were  the  solution  design 
variables  to  the  4-DV  displacement  constraint  problem  for  0.25  m.  The  use  of  another 
answer  to  begin  this  process  is  the  only  weak  point  of  this  algorithm. 

When  the  trust  region  strategy  described  in  Section  3.9  is  properly  implemented, 
the  32-DV  problem  results  shown  in  Table  4.3  are  obtained.  In  particular,  for  all  but  the 
first  case  where  the  tip  displacement  is  constrained  to  0.25  meters,  the  32-DV  results  offer 
an  average  28.7%  reduction  in  mass  compared  to  the  4-DV  results  while  satisfying  all 
constraints.  The  cost  here  is  an  average  four  times  increase  in  time  required  to  determine 
the  final  solution  over  the  4-DV  problem.  However,  these  increased  times  are  all  still  only 
between  25  and  77  seconds  total  for  the  solution.  This  level  of  design  variable  refinement 
seems  to  offer  the  best  compromise  in  the  displacement  constraint  problems  in  terms  of 
time  required,  value  of  the  objective  function,  and  the  value  of  the  limit  load,  Psafe- 

4-3  Sixty-four  Design  Variables  -  Displacement  Constraints 

Refining  the  mesh  to  a  total  of  16  beam  elements  yields  results  consistent  with 
those  from  the  previous  two  cases  as  shown  in  Tables  4.3  and  4.4  and  Figures  4.1  and 
4.2.  In  general,  the  objective  function  value  either  decreases  or  is  within  1%  of  the  value 
of  the  eight  element  mesh.  An  increase  in  computational  effort  is  present  as  the  number 
of  finite  difference  calculations  is  doubled.  The  time  cost  is  approximately  four  times 
that  of  the  eight  element  mesh  with  little  improvement  in  overall  design.  The  exception 
is  the  case  with  the  tip  displacement  constraint  set  to  0.25  meters.  It  appears  that  the 
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larger,  stiffer  models  (i.e.  higher  mass)  with  smaller  allowed  tip  displacements  may  stand 
to  lose  more  mass  with  a  more  refined  mesh. 

4-4  Four  Design  Variables  -  Stress  Constraints 

Optimization  is  performed  using  GEBT-ESL  on  an  eight  element  model  where 
the  von  Mises  stress  constraint  for  aluminum  with  a  1.5  times  safety  factor  applied 
is  179  MPa.  Cases  with  4  and  32  independent  design  variables  are  considered.  The 
least  amount  of  computational  effort  is  required  for  the  4-DV  problem,  while  the  32-DV 
problem  allowed  refinement  in  the  optimized  structure.  In  particular,  Table  4.8  shows 
for  the  32-DV  case  that  all  but  one  element  are  within  10%  of  the  applied  nonlinear 
stress  constraint.  Also  shown  is  that  elements  5  and  8  actually  exceed  the  allowed 
stress  constraint  by  5%  and  2%,  respectively.  This  discrepancy  arises  from  the  outer 
loop  convergence  criteria  used  in  the  problem.  The  convergence  tolerance  is  defined  in 
Eq.  (3.12).  This  outer  loop  tolerance  must  be  set  to  10~2  to  allow  convergence.  If  set 
to  a  tighter  tolerance,  the  optimizer  may  enter  an  endless  loop  where  the  convergence 
tolerance  e  is  never  satisfied.  With  this  relatively  “loose”  convergence  criteria,  differences 
between  the  approximate  and  actual  calculated  nonlinear  stresses  will  show  up  in  the 
second  digit.  Furthermore,  it  should  be  noted  that  all  of  the  values  presented  in  Table 
4.8  are  the  largest  von  Mises  stresses  chosen  from  values  at  the  four  corners  of  the  mid¬ 
point  cross  section.  This  method  satisfies  the  goal  of  this  research  by  showing  that  the 
constraints  can  be  applied  effectively.  With  further  refinement  of  the  approach,  tighter 
tolerances  and  evaluation  of  the  stresses  at  the  nodes  vice  the  mid  point  of  the  beam  will 
yield  even  more  accurate  designs. 

4-5  Thirty-two  Design  Variables  -  Stress  Constraints 

The  solution  to  the  32-DV  stress  constraint  problem  required  7  nonlinear  evalua¬ 
tions  and  207  seconds  to  return  the  results  shown  in  Table  4.8.  This  design  is  the  most 
flexible  one  presented  in  research  that  is  known  to  satisfy  displacement,  stress,  and  limit 
load  constraints.  Note  the  design  has  fully  stressed  members  in  elements  1  through  4 
with  a  well  behaved  taper  for  both  chord  and  height  while  skin  thickness  is  minimum 
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Table  4.8:  Optimization  results  for  beam  elements  with  stress  constraints. 


Solution 

fobj 

P safe 

^ tip 

elern 

c 

h 

ta 

tyj 

®  vm 

Method 

(kg) 

(kN) 

(m) 

(m) 

(m) 

(m) 

(m) 

(MPa) 

1 

1.0164 

0.1525 

0.0050 

0.0055 

177 

2 

1.0164 

0.1525 

0.0050 

0.0055 

121 

3 

1.0164 

0.1525 

0.0050 

0.0055 

72 

4-DV 

908.6 

127 

0.5591 

4 

1.0164 

0.1525 

0.0050 

0.0055 

33 

GEBT-ESL 

5 

1.0164 

0.1525 

0.0050 

0.0055 

32 

6 

1.0164 

0.1525 

0.0050 

0.0055 

76 

7 

1.0164 

0.1525 

0.0050 

0.0055 

131 

8 

1.0164 

0.1525 

0.0050 

0.0055 

179 

1 

1.3897 

0.2085 

0.0050 

0.0071 

179 

2 

1.1894 

0.1784 

0.0050 

0.0063 

179 

3 

0.8748 

0.1312 

0.0050 

0.0085 

180 

32-DV 

662.7 

53 

0.7771 

4 

0.5003 

0.0250 

0.0050 

0.0384 

179 

GEBT-ESL 

5 

0.5000 

0.0250 

0.0050 

0.0050 

57 

6 

0.5000 

0.0250 

0.0050 

0.0050 

90 

7 

0.5000 

0.0250 

0.0050 

0.0050 

95 

8 

0.5000 

0.0250 

0.0050 

0.0050 

140 

gage  for  all  8  elements.  The  interesting  variable  behavior  is  in  the  web  thickness  for  ele¬ 
ments  3  and  4.  The  skin  thickness  in  element  3  is  larger  that  elements  1  and  2,  although 
only  by  roughly  20  -  40%.  Element  4,  on  the  other  hand,  is  roughly  5  to  7  times  larger 
than  the  thicknesses  inboard  of  it.  From  an  aerodynamic  point  of  view,  this  is  the  easiest 
place  to  have  a  non-tapering  dimension  without  impact  to  external  planform  and  profile. 

Geometric  bend/twist  coupling  is  evident  in  Fig.  4.3.  Both  wings  have  a  discon¬ 
tinuity  at  the  3/4  span  point,  but  inboard  of  that  point  both  wings  appear  to  be  fairly 
smooth  in  their  nodal  rotations.  Considering  all  nodes  except  for  the  wing  tip  for  the  aft 
wing  shows  the  geometric  nonlinearity  that  drives  the  need  for  GEBT  as  the  nonlinear 
solver.  The  twist  in  the  aft  wing  approaches  1  degree  of  twist  at  1/4  span  before  reversing 
to  more  than  -4  degrees  of  twist  at  the  3/4  span  node.  This  change  in  sign  as  well  as  the 
magnitude  of  the  twist  indicates  the  wing  is  undergoing  considerable  deformation  apart 
from  the  calculated  tip  deflection  of  0.7771  m.  In  fact,  the  inboard  half  of  the  aft  wing 
will  be  experiencing  a  positive  increase  in  angle  of  attack  while  the  outborad  section  of 
the  aft  wing  experiences  over  4  degrees  less  angle  of  attack.  This  disparity  in  loading 
conditions  could  lead  to  flutter  or  failure  when  aerodynamic  forces  are  considered. 
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Figure  4.3:  Plot  taken  from  32- DV  stress  constraint  optimization  showing  twist  in  fore 
(blue  solid  line)  and  aft  (red  dashed  line)  wings  as  a  function  of  normalized  distance  from 
the  wing  root. 

Within  the  stress  constraint  optimization  problem  von  Mises  stresses  are  calculated 
for  32  different  points  in  an  8  element  structure.  Each  element  compares  the  maximum 
von  Mises  stress  present  to  the  stress  constraint  defined  in  Eq.  (3.10).  The  calculation  of 
the  von  Mises  stress  produces  a  vector  of  nonlinear  to  linear  stress  ratios  ot  defined  in  Eq. 
(2.10).  As  the  optimization  process  returns  to  the  outer  loop  to  update  the  nonlinear 
analysis,  a  new  ct  is  calculated.  Instead  of  converging  to  unity,  ct  does  not  converge  to 
anything  meaningful.  The  stresses  and  stress  ratio  from  the  initial  design  and  the  final 
design  for  the  32-DV  stress  constraint  problem  is  shown  in  Table  4.9. 

It  is  possible  that  the  random  nature  of  the  values  for  ot  in  Table  4.9  is  a  function 
of  using  the  displacement  based  equivalent  static  loads  and  the  stress  ratio  as  described 
by  MSC  in  their  implementation  of  ESL  [35].  It  is  unknown  at  this  time  if  using  the 
stress  equivalent  static  loads  described  by  Lee,  et  al.  will  produce  a  different  result  [5]. 
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Table  4.9:  Stress  ratio  results  for  32-DV  stress  constraint  problem 


First  Iteration 

Final  Design 

Element 

&NL-VM 

&L-VM 

OL 

<tl-vm 

&NL-VM 

gl-vm 

OL 

a  L-VM 

(MPa) 

(MPa) 

(MPa) 

(MPa) 

(MPa) 

(MPa) 

31.2 

31.2 

1.00 

31.2 

177.0 

179.0 

0.99 

177.0 

1 

21.8 

19.8 

1.10 

21.8 

105.0 

99.2 

1.06 

105.0 

1 

33.2 

32.9 

1.01 

33.2 

179.0 

175.0 

1.03 

179.0 

20.0 

21.6 

0.93 

20.0 

94.1 

94.8 

0.99 

94.1 

21.4 

22.2 

0.96 

21.4 

176.0 

216.0 

0.82 

176.0 

o 

16.6 

15.5 

1.07 

16.6 

109.0 

141.0 

0.78 

109.0 

z 

23.4 

22.1 

1.06 

23.4 

179.0 

138.0 

1.30 

179.0 

15.1 

15.4 

0.98 

15.1 

97.7 

62.8 

1.56 

97.7 

11.9 

13.4 

0.89 

11.9 

175.0 

315.0 

0.56 

175.0 

Q 

11.1 

10.9 

1.02 

11.1 

111.0 

240.0 

0.46 

111.0 

o 

13.9 

11.6 

1.20 

13.9 

180.0 

39.6 

4.54 

180.0 

10.0 

9.1 

1.11 

10.0 

96.7 

37.5 

2.58 

96.7 

3.8 

4.2 

0.90 

3.8 

179.0 

774.0 

0.23 

179.0 

A 

5.7 

5.5 

1.03 

5.7 

150.0 

729.0 

0.21 

150.0 

4 

5.3 

2.0 

2.73 

5.3 

174.0 

382.0 

0.45 

174.0 

5.7 

3.2 

1.80 

5.7 

134.0 

427.0 

0.31 

134.0 

5.7 

7.8 

0.73 

5.7 

20.6 

189.0 

0.11 

20.6 

5 

3.6 

6.2 

0.58 

3.6 

38.4 

243.0 

0.16 

38.4 

5.7 

1.9 

3.00 

5.7 

44.1 

170.0 

0.26 

44.1 

4.9 

1.7 

2.85 

4.9 

56.8 

115.0 

0.49 

56.8 

10.9 

12.9 

0.85 

10.9 

89.5 

138.0 

0.65 

89.5 

6 

12.1 

15.8 

0.76 

12.1 

18.4 

104.0 

0.18 

18.4 

9.9 

6.8 

1.45 

9.9 

49.8 

50.2 

0.99 

49.8 

14.1 

9.7 

1.45 

14.1 

60.0 

83.6 

0.72 

60.0 

16.3 

17.4 

0.94 

16.3 

94.5 

330.0 

0.29 

94.5 

7 

21.6 

24.6 

0.88 

21.6 

57.5 

318.0 

0.18 

57.5 

15.0 

13.2 

1.14 

15.0 

68.6 

179.0 

0.38 

68.6 

23.8 

20.4 

1.17 

23.8 

79.9 

191.0 

0.42 

79.9 

22.0 

22.2 

0.99 

22.0 

78.3 

152.0 

0.52 

78.3 

O 

31.0 

33.1 

0.94 

31.0 

106.0 

206.0 

0.51 

106.0 

O 

20.6 

19.8 

1.04 

20.6 

55.4 

23.5 

2.35 

55.4 

33.3 

30.8 

1.08 

33.3 

140.0 

30.7 

4.56 

140.0 
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4-6  Failure  to  Converge  to  Displacement  Constraint 

There  exist  some  combinations  of  move  limits,  move  reductions,  and  allowable  tip 
deflection  that  produce  results  that  do  not  converge  according  to  Eq.  (3.12).  An  example 
of  this  is  shown  in  Fig.  4.4.  Due  to  the  nonlinearity  of  the  joined-wing  problem,  it  is 
possible  for  the  linear  solution  to  converge  within  the  inner  loop  after  making  changes  in 
design  variables  that  have  greater  sensitivity  than  the  linear  anlaysis  is  able  to  capture, 
even  if  a  trust  region  strategy  is  employed.  As  a  result,  the  ensuing  outer  loop  nonlinear 
analysis  with  GEBT  provides  a  tip  displacement  that  is  significantly  different  than  the 
value  at  the  end  of  the  inner  loop  optimization.  The  disconnect  then  comes  with  the 
reduction  in  move  limits  before  the  next  inner  loop  begins.  When  the  step  size  is  too 
small,  the  inner  loop  will  not  be  able  to  move  away  from  the  “bad”  set  of  design  variables 
enough  to  get  out  of  the  region  in  the  design  space  where  the  nonlinear  analysis  is  so 
sensitive.  A  possible  solution  could  be  to  apply  the  displacement  constraint  as  a  specific 
criterion  in  the  trust  region  strategy  in  conjunction  with  the  merit  function  shown  in 
Eqs.  (3.21)  and  (3.22). 

4-7  Convergence  from  Outside  the  Limit  Load  Feasible  Design  Space 

After  applying  the  limit  load  approximation  technique  described  in  Section  3.8.2, 
a  test  was  designed  to  determine  how  well  the  approximation  technique  dealt  with  a 
design  that  fails  the  limit  load  test.  The  chosen  design  was  the  lower  limit,  or  minimum 
gage,  for  all  design  variables  as  shown  in  Eq.  (3.10)  in  a  32-DV  displacement  constraint 
problem.  The  allowable  tip  deflection,  uzauow,  was  0.70  m. 

The  first  issue  that  arose  was  the  trust  region  strategy  described  in  Section  3.9 
worked  as  desired.  This  meant  that  even  though  the  linear  response  optimization  at¬ 
tempted  to  increase  the  size  of  the  design  variables  to  produce  a  more  robust  structure, 
the  trust  region  strategy  negated  these  attempts  due  to  a  failure  on  the  part  of  the  new 
solution  to  improve  upon  the  mass  of  the  structure.  Since  the  design  is  at  a  minimum 
possible  mass  of  398.94  kg  when  all  design  variables  are  minimum  gage,  the  TR  strategy 
believes  it  has  the  best  solution  available,  regardless  of  the  optmimzers  efforts  to  achieve 
a  feasible  design. 
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Figure  4.4:  Three  plots  showing  behavior  of  32-DV  nonlinear  wing  tip  deflection  at¬ 
tempting  to  converge  to  uup  =  1.15  m,  mass,  and  value  of  the  the  tolerance,  Equation 
(3.12) 


With  the  trust  region  strategy  turned  off,  the  optimizer  was  able  to  change  the 
design  variables  in  a  manner  that  increased  the  stiffness  and  mass  of  the  design.  However, 
as  stated  previously,  the  convergence  of  this  problem  to  a  feasible  design  is  still  sensitive 
to  the  chosen  move  limits  and  move  reduction  ratio  defined  in  the  problem.  Move  limits 
of  0.5  and  1.5  applied  to  all  design  variables  were  too  small  and  too  large,  respectively 
when  paried  with  a  70%  move  reduction.  However,  when  a  move  limit  of  1.0  is  combined 
with  a  50%  move  reduction,  the  design  shown  in  Table  4.10  is  the  result.  This  design 
is  significant  in  that  fQbj  =  726.1  kg,  and  Psafe  =  67  kN,  a  feasible  design.  It  does  not, 
however,  lie  on  the  Pareto  curve  for  the  32-DV  displacement  problem  shown  in  Fig.  4.2. 
The  estimated  value  of  fobj  is  652  kg  for  a  displacement  constraint  problem  where  uzauow 
=  0.70  m.  The  other  interesting  characteristic  is  that  the  optimizer  chose  a  path  that 
results  in  both  wings  tapering,  instead  of  one  going  to  minimum  gage  as  is  common 
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Tabic  4.10:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uauow  =  0.70 
m.  Limit  load  estimated  using  ratio  method  described  in  Section  3.8.2. 


Element 

c 

(m) 

h 

(m) 

^ skin 

(m) 

^ web 

(m) 

1 

1.0263 

0.1539 

0.0050 

0.0050 

2 

0.8759 

0.1211 

0.0050 

0.0050 

3 

0.6653 

0.0879 

0.0050 

0.0050 

4 

0.5000 

0.0607 

0.0050 

0.0050 

5 

0.5005 

0.0751 

0.0050 

0.0050 

6 

0.7495 

0.1124 

0.0050 

0.0050 

7 

1.0512 

0.1577 

0.0050 

0.0050 

8 

1.2630 

0.1783 

0.0050 

0.0050 

when  starting  from  the  initial  design  variable  values  shown  in  Table  3.1.  It  is  possible 
that  more  experimentation  with  move  limits  and  move  reduction  ratio  would  produce  an 
optimized  design  that  lies  on  the  Pareto  front  while  remaining  feasible  with  respect  to 
the  limit  load  and  starting  from  an  infeasible  design.  However,  the  larger  issue  lies  with 
the  incompatibillity  of  the  trust  region  with  the  minimum  gage  initial  design.  Either 
a  means  of  determining  when  to  use  the  trust  region  needs  to  be  defined,  or  all  initial 
design  variables  should  be  chosen  such  that  the  initial  design  has  greater  mass  than  the 
anticipated  final  design. 
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V.  Conclusions 


5.1  Overview  of  Research  Effort 

A  full  nonlinear  response  optimization  on  the  joined-wing  design  is  presented.  De¬ 
sign  variables  relating  to  the  cross  section  of  a  simple  box  beam  were  sized  using 
an  SQP  optimization  routine.  Geometrically  exact  beam  theory  is  shown  to  be  an  effec¬ 
tive  tool  in  modeling  the  geometric  nonlinearities  exhibited  by  the  joined-wing  design. 
Further  discretization  of  the  fore  and  aft  wings  to  a  total  of  16  elements  with  64  de¬ 
sign  variables  did  not  significantly  improve  on  the  optimized  design  in  the  displacement 
constraint  problem.  Considering  this  and  the  fact  that  the  time  required  increased  on 
average  four  times,  it  is  clear  that  the  8  element,  32  design  variable  model  offers  the  best 
compromise  in  time  required  and  final  optimized  mass  of  the  structure. 

5.2  Conclusions 

The  use  of  equivalent  static  loads  with  GEBT  are  shown  to  be  an  effective  solution 
method  while  offering  significant  savings  in  computational  effort  required.  There  exist 
unexplained  sensitivities  in  the  use  of  displacement  based  equivalent  static  loads  that 
occasionally  permit  the  optimization  routine  to  converge  to  a  displacement  constraint 
in  the  nonlinear  solution  that  is  not  the  same  as  the  applied  constraint  in  the  linear 
approximation,  ft  is  believed  that  these  sensitivities  are  due  at  least  in  part  to  the 
move  limit  and  move  reduction  allowed  from  iteration  to  iteration  when  a  trust  region 
strategy  is  applied.  The  current  research  has  only  applied  the  trust  region  strategy  to  the 
displacement  constraint  problem,  so  evaluation  of  the  trust  region  effects  do  not  apply 
to  stress  constraint  problems  at  this  time. 

The  need  for  consideration  of  additional  constraints  is  shown.  In  particular,  the 
limit  load  criterion  and  the  requirement  that  the  design  remain  feasible  with  respect  to 
the  limit  load  is  a  key  safety  consideration.  While  a  finite  difference  approach  to  sensitiv¬ 
ity  analysis  is  used,  the  determination  of  this  constraint  is  problematic.  The  easiest  and 
computationally  costliest  approach  is  to  recompute  the  load-deflection  curve  whenever  a 
change  is  made  to  the  design  variables.  This  method  renders  the  use  of  equivalent  static 
loads  ineffective  due  to  the  number  of  nonlinear  calculations  required  to  determine  the 
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load-deflection  curve  for  each  iteration  in  the  linear  response  optimization.  An  alternate 
means  of  determining  the  limit  load  using  a  second  order  polynomial  approximation 
which  updates  the  slope  term  during  the  linear  FEM  optimization  was  presented  and 
shown  to  be  effective. 

5.3  Future  Work 

Future  work  should  focus  on  determining  the  how  the  move  limit  and  move  reduc¬ 
tion  sizes  affect  the  solution  so  that  more  general  optimization  problems  can  be  set  up 
with  less  interaction  required  on  the  part  of  the  user.  This  will  require  a  more  detailed 
analysis  of  the  logic  employed  in  the  trust  region  strategy  currently  employed.  It  is  also 
planned  that  the  trust  ratio  would  be  expanded  to  include  stress  constraint  problems  as 
well  as  displacement  constraints. 

Also,  working  to  implement  the  stress  based  equivalent  static  loads  (SBESL)  de¬ 
scribed  by  Lee  et  al.  [5]  and  making  a  comparison  between  SBESL  and  the  results 
presented  in  this  research  would  be  beneficial.  While  MSC  has  opted  for  the  simpler 
implementation  shown  here,  it  remains  to  be  seen  if  there  is  an  appreciable  difference 
in  results  when  using  one  method  or  the  other.  The  expectation  at  this  time  is  that  a 
would  reduce  to  unity  when  the  optimizatino  returns  a  final  solution.  If  there  is  little  to 
no  difference  between  SBESL  and  the  method  used  in  this  research,  it  begs  the  question 
of  how  important  is  the  source  data  for  the  stress  case  as  long  as  a  stress  ratio  correction 
is  applied. 

A  scale  model  of  an  optimized  design  made  of  aluminum  box  tubes  is  planned  for 
laboratory  experimentation  to  compare  and  contrast  to  the  results  obtained  using  the 
method  outlined  in  this  research.  Tip  loads  and  joint  loads  are  planned  on  a  design  that 
locates  the  wing  joint  at  less  than  100%  span  of  the  forward  wing. 

Finally,  incorporation  of  shape  optimization  of  planform  will  allow  a  greater  level 
of  optimization  based  on  aerodynamic,  structural,  and  sensor  effectiveness  sensitivities. 
Consideration  for  this  implementation  must  include  the  application  of  coordinate  trans¬ 
formation  from  a  change  in  angles  in  both  sweep  and  dihedral  for  both  wings.  With  these 
design  variables  included,  a  very  complete  and  robust  design  analysis  can  be  conducted 
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on  airfoil  shapes  and  orientation  which  should  show  the  path  to  the  most  efficient  design 
for  the  joined-wing  Sensorcraft  concept. 
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Appendix  A.  Optimization  Results  for  JrDV  Optimization  with 

Displacement  Constraints 


THE  following  tables  show  the  cross  sectional  dimensions  of  the  final  design  obtained 
for  4,  32,  and  64-DV  problems  with  varying  displacement  constraints  applied. 


Table  A.l:  Design  variables  by  element  for  4-DV  GEBT-ESL  design,  P- limit  estimation 
used  for  limit  load  calculation. 


fallow 

c 

(m) 

h 

(m) 

t'skin 

(m) 

t-web 

(m) 

0.25 

1.3368 

0.2005 

0.0050 

0.0050 

0.35 

1.1954 

0.1793 

0.0050 

0.0050 

0.45 

1.0010 

0.1650 

0.0050 

0.0050 

0.55 

1.0299 

0.1545 

0.0050 

0.0050 

0.65 

0.9753 

0.1463 

0.0050 

0.0050 

0.75 

0.9311 

0.1397 

0.0050 

0.0050 

0.85 

0.8945 

0.1342 

0.0050 

0.0050 

0.95 

0.8634 

0.1295 

0.0050 

0.0050 

1.05 

0.8367 

0.1255 

0.0050 

0.0050 

1.15 

0.8135 

0.1220 

0.0050 

0.0050 

1.25 

0.7931 

0.1190 

0.0050 

0.0050 
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Appendix  B.  Optimization  Results  for  32-DV  Optimization  with 

Displacement  Constraints 


THE  following  tables  show  the  cross  sectional  dimensions  of  the  final  design  obtained 
for  4,  32,  and  64-DV  problems  with  varying  displacement  constraints  applied. 


Table  B.l:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uaiiow  =  0.25 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

t web 

(m) 

1 

1.3969 

0.1777 

0.0050 

0.0050 

2 

1.3969 

0.1069 

0.0050 

0.0050 

3 

1.3969 

0.0719 

0.0050 

0.0050 

4 

1.3969 

0.0719 

0.0050 

0.0050 

5 

1.3969 

0.0719 

0.0050 

0.0050 

6 

1.3969 

0.1063 

0.0050 

0.0050 

7 

1.3969 

0.1759 

0.0050 

0.0050 

8 

1.4818 

0.2218 

0.0050 

0.0050 

Table  B.2:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uaiiow  =  0.35 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

tweb 

(m) 

1 

0.5775 

0.0300 

0.0050 

0.0050 

2 

0.5775 

0.0300 

0.0050 

0.0050 

3 

0.5775 

0.0300 

0.0050 

0.0050 

4 

0.5775 

0.0300 

0.0050 

0.0050 

5 

0.7161 

0.1074 

0.0050 

0.0050 

6 

1.1824 

0.1770 

0.0050 

0.0050 

7 

1.5100 

0.2260 

0.0050 

0.0050 

8 

1.7775 

0.2660 

0.0050 

0.0050 
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Table  B.3:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uauow  =  0.45 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

t skin 

(m) 

t web 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.6587 

0.0985 

0.0050 

0.0050 

6 

1.0924 

0.1639 

0.0050 

0.0050 

7 

1.3913 

0.2087 

0.0050 

0.0050 

8 

1.6363 

0.2449 

0.0050 

0.0050 

Table  B.4:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uaiiow  =  0.55 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

tweb 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.6192 

0.0922 

0.0050 

0.0050 

6 

1.0186 

0.1525 

0.0050 

0.0050 

7 

1.2963 

0.1942 

0.0050 

0.0050 

8 

1.5238 

0.2274 

0.0050 

0.0050 

Table  B.5:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uaiiow  =  0.65 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

tweb 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5939 

0.0888 

0.0050 

0.0050 

6 

0.9630 

0.1442 

0.0050 

0.0050 

7 

1.2206 

0.1828 

0.0050 

0.0050 

8 

1.4313 

0.2145 

0.0050 

0.0051 
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Table  B.6:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uauow  =  0.75 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

t skin 

(m) 

t web 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5634 

0.0845 

0.0050 

0.0050 

6 

0.9197 

0.1380 

0.0050 

0.0050 

7 

1.1670 

0.1750 

0.0050 

0.0050 

8 

1.3583 

0.2037 

0.0050 

0.0050 

Table  B.7:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uaiiow  =  0.85 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

tweb 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5440 

0.0816 

0.0050 

0.0050 

6 

0.8815 

0.1322 

0.0050 

0.0050 

7 

1.1174 

0.1676 

0.0050 

0.0050 

8 

1.3011 

0.1952 

0.0050 

0.0050 

Table  B.8:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uaiiow  =  0.95 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

tweb 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5355 

0.0803 

0.0050 

0.0050 

6 

0.8484 

0.1273 

0.0050 

0.0050 

7 

1.0696 

0.1604 

0.0050 

0.0050 

8 

1.2549 

0.1882 

0.0050 

0.0050 
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Table  B.9:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uauow  =  1.05 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

t skin 

(m) 

t web 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5163 

0.0774 

0.0050 

0.0050 

6 

0.8301 

0.1245 

0.0050 

0.0050 

7 

1.0386 

0.1558 

0.0050 

0.0050 

8 

1.2033 

0.1805 

0.0050 

0.0050 

Table  B.10:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uauow  =  1.15 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

t web 

(m) 

1 

0.5000 

0.0265 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5119 

0.0762 

0.0050 

0.0050 

6 

0.8113 

0.1217 

0.0050 

0.0050 

7 

1.0017 

0.1500 

0.0050 

0.0050 

8 

1.1669 

0.1731 

0.0050 

0.0050 

Table  B.ll:  Design  variables  by  element  for  32-DV  GEBT-ESL  design,  uaiiow  =  1-25 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

tweb 

(m) 

1 

0.8196 

0.1220 

0.0050 

0.0050 

2 

0.6969 

0.0957 

0.0050 

0.0050 

3 

0.6969 

0.0394 

0.0050 

0.0050 

4 

0.6969 

0.0349 

0.0050 

0.0050 

5 

0.6969 

0.0375 

0.0050 

0.0050 

6 

0.7245 

0.1087 

0.0050 

0.0050 

7 

0.8531 

0.1231 

0.0050 

0.0050 

8 

0.9362 

0.1356 

0.0050 

0.0050 
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Appendix  C.  Optimization  Results  for  64-DV  Optimization  with 

Displacement  Constraints 


THE  following  tables  show  the  cross  sectional  dimensions  of  the  final  design  obtained 
for  64-DV  problems  with  varying  displacement  constraints  applied. 


Table  C.l:  Design  variables  by  element  for  64-DV  GEBT-ESL  design,  uaiiow  =  0.25 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

^ web 

(m) 

1 

0.5000 

0.0288 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5000 

0.0250 

0.0050 

0.0050 

6 

0.5000 

0.0250 

0.0050 

0.0050 

7 

0.5000 

0.0250 

0.0050 

0.0050 

8 

0.5000 

0.0250 

0.0050 

0.0050 

9 

0.5902 

0.0885 

0.0050 

0.0050 

10 

0.9638 

0.1446 

0.0050 

0.0050 

11 

1.2166 

0.1825 

0.0050 

0.0050 

12 

1.4304 

0.2146 

0.0050 

0.0050 

13 

1.6155 

0.2423 

0.0050 

0.0050 

14 

1.7868 

0.2680 

0.0050 

0.0050 

15 

1.9419 

0.2910 

0.0050 

0.0050 

16 

2.0895 

0.3134 

0.0050 

0.0050 
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Tabic  C.2:  Design  variables  by  element  for  64-DV  GEBT-ESL  design,  uauow  =  0.35 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

t'skin 

(m) 

t"web 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5000 

0.0250 

0.0050 

0.0050 

6 

0.5000 

0.0250 

0.0050 

0.0050 

7 

0.5000 

0.0250 

0.0050 

0.0050 

8 

0.5000 

0.0250 

0.0050 

0.0050 

9 

0.5198 

0.0780 

0.0050 

0.0050 

10 

0.8544 

0.1282 

0.0050 

0.0050 

11 

1.0881 

0.1632 

0.0050 

0.0050 

12 

1.2801 

0.1920 

0.0050 

0.0050 

13 

1.4434 

0.2165 

0.0050 

0.0050 

14 

1.5896 

0.2384 

0.0050 

0.0050 

15 

1.7251 

0.2588 

0.0050 

0.0050 

16 

1.8452 

0.2768 

0.0050 

0.0050 

Table  C.3:  Design  variables  by  element  for  64-DV  GEBT-ESL  design,  uauow  =  0.45 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

t web 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5000 

0.0250 

0.0050 

0.0050 

6 

0.5000 

0.0250 

0.0050 

0.0050 

7 

0.5000 

0.0250 

0.0050 

0.0050 

8 

0.5000 

0.0250 

0.0050 

0.0050 

9 

0.5000 

0.0750 

0.0050 

0.0050 

10 

0.7840 

0.1176 

0.0050 

0.0050 

11 

0.9975 

0.1496 

0.0050 

0.0050 

12 

1.1669 

0.1750 

0.0050 

0.0050 

13 

1.3196 

0.1979 

0.0050 

0.0050 

14 

1.4535 

0.2180 

0.0050 

0.0050 

15 

1.5741 

0.2361 

0.0050 

0.0050 

16 

1.6860 

0.2529 

0.0050 

0.0050 

60 


Table  C.4:  Design  variables  by  element  for  64-DV  GEBT-ESL  design,  uauow  =  0.55 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

t"web 

(m) 

1 

1.6297 

0.2442 

0.0050 

0.0050 

2 

1.5151 

0.2273 

0.0050 

0.0050 

3 

1.3929 

0.2089 

0.0050 

0.0050 

4 

1.2607 

0.1891 

0.0050 

0.0050 

5 

1.1120 

0.1668 

0.0050 

0.0050 

6 

0.9427 

0.1411 

0.0050 

0.0050 

7 

0.7385 

0.1096 

0.0050 

0.0050 

8 

0.5000 

0.0486 

0.0050 

0.0050 

9 

0.5000 

0.0250 

0.0050 

0.0050 

10 

0.5000 

0.0250 

0.0050 

0.0050 

11 

0.5000 

0.0250 

0.0050 

0.0050 

12 

0.5000 

0.0250 

0.0050 

0.0050 

13 

0.5000 

0.0250 

0.0050 

0.0050 

14 

0.5000 

0.0250 

0.0050 

0.0050 

15 

0.5000 

0.0250 

0.0050 

0.0050 

16 

0.5000 

0.0250 

0.0050 

0.0050 

Table  C.5:  Design  variables  by  element  for  64-DV  GEBT-ESL  design,  uauow  =  0.65 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

t web 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5000 

0.0250 

0.0050 

0.0050 

6 

0.5000 

0.0250 

0.0050 

0.0050 

7 

0.5000 

0.0250 

0.0050 

0.0050 

8 

0.5000 

0.0250 

0.0050 

0.0050 

9 

0.5000 

0.0597 

0.0050 

0.0050 

10 

0.7009 

0.1046 

0.0050 

0.0050 

11 

0.8825 

0.1324 

0.0050 

0.0050 

12 

1.0311 

0.1547 

0.0050 

0.0050 

13 

1.1613 

0.1742 

0.0050 

0.0050 

14 

1.2751 

0.1912 

0.0050 

0.0050 

15 

1.3814 

0.2065 

0.0050 

0.0050 

16 

1.4800 

0.2205 

0.0050 

0.0050 
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Table  C.6:  Design  variables  by  element  for  64-DV  GEBT-ESL  design,  uauow  =  0.75 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

t'skin 

(m) 

t"web 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5000 

0.0250 

0.0050 

0.0050 

6 

0.5000 

0.0250 

0.0050 

0.0050 

7 

0.5000 

0.0250 

0.0050 

0.0050 

8 

0.5000 

0.0250 

0.0050 

0.0050 

9 

0.5000 

0.0540 

0.0050 

0.0050 

10 

0.6756 

0.1001 

0.0050 

0.0050 

11 

0.8425 

0.1264 

0.0050 

0.0050 

12 

0.9842 

0.1476 

0.0050 

0.0050 

13 

1.1056 

0.1658 

0.0050 

0.0050 

14 

1.2123 

0.1819 

0.0050 

0.0050 

15 

1.3162 

0.1963 

0.0050 

0.0050 

16 

1.4012 

0.2094 

0.0050 

0.0050 

Table  C.7:  Design  variables  by  element  for  64-DV  GEBT-ESL  design,  uauow  =  0.85 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(nr) 

h 

(nr) 

tskin 

(nr) 

t web 

(nr) 

1 

0.5001 

0.0250 

0.0050 

0.0050 

2 

0.5001 

0.0250 

0.0050 

0.0050 

3 

0.5001 

0.0250 

0.0050 

0.0050 

4 

0.5001 

0.0250 

0.0050 

0.0050 

5 

0.5001 

0.0250 

0.0050 

0.0050 

6 

0.5001 

0.0250 

0.0050 

0.0050 

7 

0.5000 

0.0250 

0.0050 

0.0050 

8 

0.5001 

0.0250 

0.0050 

0.0050 

9 

0.5001 

0.0494 

0.0050 

0.0050 

10 

0.6511 

0.0970 

0.0050 

0.0050 

11 

0.8134 

0.1220 

0.0050 

0.0050 

12 

0.9428 

0.1414 

0.0050 

0.0050 

13 

1.0627 

0.1594 

0.0050 

0.0050 

14 

1.1627 

0.1744 

0.0050 

0.0050 

15 

1.2533 

0.1878 

0.0050 

0.0050 

16 

1.3425 

0.1999 

0.0050 

0.0050 
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Table  C.8:  Design  variables  by  element  for  64-DV  GEBT-ESL  design,  uauow  =  0.95 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

t"web 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5000 

0.0250 

0.0050 

0.0050 

6 

0.5000 

0.0250 

0.0050 

0.0050 

7 

0.5000 

0.0250 

0.0050 

0.0050 

8 

0.5000 

0.0250 

0.0050 

0.0050 

9 

0.5000 

0.0461 

0.0050 

0.0050 

10 

0.6349 

0.0941 

0.0050 

0.0050 

11 

0.7871 

0.1179 

0.0050 

0.0050 

12 

0.9142 

0.1371 

0.0050 

0.0050 

13 

1.0235 

0.1535 

0.0050 

0.0050 

14 

1.1191 

0.1679 

0.0050 

0.0050 

15 

1.2077 

0.1806 

0.0050 

0.0050 

16 

1.2923 

0.1921 

0.0050 

0.0050 

Table  C.9:  Design  variables  by  element  for  64-DV  GEBT-ESL  design,  uauow  =  1.05 
m.  No  move  limits  or  limit  load  calculation. 


Element 

c 

(m) 

h 

(m) 

tskin 

(m) 

t web 

(m) 

1 

0.5000 

0.0250 

0.0050 

0.0050 

2 

0.5000 

0.0250 

0.0050 

0.0050 

3 

0.5000 

0.0250 

0.0050 

0.0050 

4 

0.5000 

0.0250 

0.0050 

0.0050 

5 

0.5000 

0.0250 

0.0050 

0.0050 

6 

0.5000 

0.0250 

0.0050 

0.0050 

7 

0.5000 

0.0250 

0.0050 

0.0050 

8 

0.5000 

0.0250 

0.0050 

0.0050 

9 

0.5000 

0.0424 

0.0050 

0.0050 

10 

0.6225 

0.0934 

0.0050 

0.0050 

11 

0.7738 

0.1161 

0.0050 

0.0050 

12 

0.8820 

0.1323 

0.0050 

0.0050 

13 

0.9966 

0.1495 

0.0050 

0.0050 

14 

1.0873 

0.1631 

0.0050 

0.0050 

15 

1.1657 

0.1749 

0.0050 

0.0050 

16 

1.2363 

0.1854 

0.0050 

0.0050 
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Appendix  D.  Third  Possible  Limit  Load  Estimation  Technique 


The  discussion  on  limit  loads  and  ways  to  estimate  them  in  Sections  3.7  and  3.8.2 
covered  approaches  used  in  this  research.  A  potential  third  option  exists  which 
is  based  on  the  sensitivity  of  the  limit  loads  with  respect  to  the  design  variables.  This  is 
shown  as 

Plim  =  Plim  +  V  P^mAX  (D.l) 

where  Pnrn  is  the  approximate  limit  load,  Pjtm  is  the  limit  load  about  which  the  sensitiv¬ 
ities  exist,  T P um  is  the  vector  of  partial  derivatives  of  Pum  with  respect  to  each  design 
variable,  and  AA"  is  the  vector  of  change  in  the  value  of  the  design  variables.  Analytic 
sensitivities  of  the  limit  load  could  be  easily  used  with  the  other  information  already 
available  to  provide  a  quick  estimate  of  the  limit  load.  Calculation  of  analytic  sensitivi¬ 
ties  should  be  included  in  future  work  as  a  tool  to  further  refine  the  optimization  of  the 
joined-wing  beam  model  used  in  this  research. 


64 


Bibliography 

1.  Bowers,  P.  M.,  Unconventional  Aircraft,  Tab  Books,  Inc.,  Blue  Ridge  Summit,  PA, 
1984. 

2.  Wolkovitch,  J.,  “Joined  Wing  Aircraft,”  U.S.  Patent  4,365,773,  December  28,  1982. 

3.  Wolkovitch,  J.,  “The  Joined-Wing:  An  Overview,”  Journal  of  Aircraft,  Vol.  23, 
No.  3,  1986,  pp.  161-178. 

4.  Gallman,  J.  W.  and  Kroo,  I.  M.  ,  “Structural  Optimization  of  Joined-Wing  Synthe¬ 
sis,”  Journal  of  Aircraft,  Vol.  33,  No.  1,  1996,  pp.  214-223. 

5.  Lee,  H.,  Kim,  Y.  I.,  and  Park  G.  J.,  et  al.,  “Structural  Optimization  of  a  Joined- 
Wing  Using  Equivalent  Static  Loads,”  Journal  of  Aircraft,  Vol.  44,  No.  4,  2007, 
pp.  1302-1308. 

6.  Kim,  Y.,  Park,  G.  J.,  Kolonay,  R.,  Blair,  M.,  and  Canfield,  R.  A.,  “Nonlinear  Re¬ 
sponse  Structural  Optimization  of  a  Joined-Wing  Using  Equivalent  Static  Loads,” 
AIAA  Journal,  Vol.  46,  No.  11,  2008,  pp.  2703-2713. 

7.  Hodges,  D.  H.,  Atilgan,  A.,  Cesnik,  C.  E.  S.,  and  Fulton,  M.,  “On  a  Simplified 
Strain  Energy  Function  for  Geometrically  Nonlinear  Behavior  of  Aniostropic  Beams,” 
Composites  Engineering,  Vol.  2,  1992,  pp.  513-526. 

8.  Cesnik,  C.  E.  S.  and  Hodges,  D.  H.,  “VABS:  A  New  Concept  for  Composite  Ro¬ 
tor  Blade  Cross  Sectional  Modeling,”  Journal  of  the  American  Helicopter  Society, 
Vol.  42,  1997,  pp.  27-38. 

9.  Yu,  W.,  Hodges,  D.  H.,  Volovoi,  V.  V.,  and  Cesnik,  C.  E.  S.,  “On  Timoshenko-like 
Modeling  of  Initially  Curved  and  Twisted  Composite  Beams,”  International  Journal 
of  Solids  and  Structures,  Vol.  39,  2002,  pp.  5101-5121. 

10.  Blair,  M.  and  Stritz,  A.,  “Finite  Element  Beam  Assemblies  with  Geometric  Bend- 
Twist  Coupling,”  AIAA/ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics, 
and  Materials  Conference,  Schaumberg,  IL,  2008. 

11.  Hodges,  D.  H.,  “A  Mixed  Variational  Formulation  Based  on  Exact  Intrinsic  Equa¬ 
tions  for  Dynamics  of  Moving  Beams,”  International  Journal  of  Solids  and  Struc¬ 
tures,  Vol.  26,  1990,  pp.  1253-1273. 

12.  Hodges,  D.  H.,  Nonlinear  Composite  Beam  Theory ,  American  Institute  of  Aeronau¬ 
tics  and  Astronautics,  Inc.,  Reston,  VA,  2006. 

13.  Kunz,  D.,  “Survey  and  Comparison  of  Engineering  Beam  Theories  for  Helicopter 
Rotor  Blades,”  Journal  of  Aircraft,  Vol.  31,  1994,  pp.  473-479. 

14.  Li,  L.,  Volovoi,  V.  V.,  and  Hodges,  D.  H.,  “Cross  Setional  Design  of  Composite  Rotor 
Blades,”  Journal  of  the  American  Helicopter  Society,  Vol.  53,  2008,  pp.  240-251. 

15.  Volovoi,  V.  V.,  Hodges,  D.  H.,  Cesnik,  C.  E.  S.,  and  Popescu,  B.,  “Assessment  of 
Beam  Modeling  Methods  for  Rotor  Blade  Applications,”  Mathematical  and  Com¬ 
puter  Modeling,  Vol.  33,  2001,  pp.  1099-1112. 


65 


16.  Johnson,  W.,  “Rotorcraft  Aerodynamics  Models  for  a  Comprehensive  Analysis,”  Pro¬ 
ceedings  of  the  54th  Annual  Forum  of  the  American  Helicopter  Society,  Washington, 
D.C.,  1998. 

17.  Johnson,  W.,  “Rotorcraft  Dynamics  Models  for  a  Comprehensive  Analysis,”  Pro¬ 
ceedings  of  the  54th  Annual  Forum  of  the  American  Helicopter  Society,  Washington, 
D.C.,  1998. 

18.  Bauchau,  O.,  “Computational  Schemes  for  Flexible,  Nonlinear  Multi- Body  Systems,” 
Multibody  System  Dynamics ,  Vol.  2,  1998,  pp.  169-225. 

19.  Saberi,  H.  A.,  Khoshlahjeh,  M.,  Ormiston,  R.  A.,  and  Rutkowski,  M.  J.,  “RCAS 
Overview  and  Application  to  Advanced  Rotorcraft  Problems,”  4th  Decennial  Spe¬ 
cialists’  Conference  on  Aeromechanics,  American  Helicopter  Society,  21-23  Jan  2004. 

20.  Cesnik,  C.  E.  S.  and  Brown,  E.  L.,  “Active  Warping  Control  of  a  Joined  Wing  Air¬ 
plane  Configuration,”  AIAA/ASME/ASCE/AHS/ASC  Structures,  Structural  Dy¬ 
namics,  and  Materials  Conference,  Hampton,  VA,  2003. 

21.  Cesnik,  C.  E.  S.  and  Ortega- Morales,  M.,  “Active  Aeroelastic  Tailoring  of  Slender 
Flexible  Wings,”  International  Forum  on  Aeroelasticity  and  Structural  Dynamics, 
Madrid,  Spain,  2001. 

22.  Blair,  M.  and  Canfield,  R.,  “A  Joined- Wing  Structural  Weight  Modeling  Study,” 
AIAA/ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics,  and  Materials 
Conference,  Denver,  CO,  Apr  2002. 

23.  Danielson,  D.  A.  and  Hodges,  D.  H.,  “Nonlinear  Beam  Kinematics  by  Decomposition 
of  the  Rotation  Tensor,”  Journal  of  Applied  Mechanics ,  Vol.  54,  1987,  pp.  258-262. 

24.  Danielson,  D.  A.  and  Hodges,  D.  H.,  “A  Beam  Theory  for  Large  Global  Rotation, 
Moderate  Local  Rotation,  and  Small  Strain,”  Journal  of  Applied  Mechanics ,  Vol.  55, 
1988,  pp.  179-184. 

25.  Parker,  D.  F.,  “The  Role  of  Saint  Venant’s  Solutions  in  Rod  and  Beam  Theories,” 
Journal  of  Applied  Mechanics ,  Vol.  46,  1979,  pp.  861-866. 

26.  Atilgan,  A.  R.  and  Hodges,  D.  H.,  “A  geometrically  nonlinear  analysis  for  nonho- 
mogeneous,  anisotropic  beams,”  Proc.  30th  AIAA/ASME/ASCE/AHS/ASC  Struc¬ 
tures,  Structural  Dynamics,  and  Materials  Conference,  Mobile,  AL,  3-5  April  1989. 

27.  Reissner,  E.,  “On  one-dimensional  large-displacement  finite-strain  beam  theory,” 
Studies  in  Applied  Mathematics ,  Vol.  52,  1973,  pp.  87-95. 

28.  Reissner,  E.,  “On  finite  deformations  of  space-curved  beams,”  ZAMP ,  Vol.  32,  1981, 
pp.  734-744. 

29.  Borri,  M.  and  Mantegazza,  P.,  “Some  Contributions  on  Structural  and  Dynamic 
Modeling  of  Helicopter  Rotor  Blades,”  lAerotecnica  Missili  e  Spazio,  Vol.  64,  1985, 
pp.  143-154. 

30.  Bachau,  O.  and  Kang,  N.,  “A  Multibody  Formulation  for  Helicopter  Structural  Dy¬ 
namic  Analysis,”  Journal  of  the  American  Helicopter  Society ,  Vol.  38,  1993,  pp.  3-14. 

31.  Berdichevskii,  V.  L.,  “Equations  of  the  theory  of  anisotropic  inhomogeneous  rods,” 
Dokl.  Akad.  Nauk.  SSR,  Vol.  228,  1976,  pp.  558-561. 


66 


32.  Popescu,  B.  and  Hodges,  D.  H.,  “On  asymptotically  correct  Timoshenko-like 
anisotropic  beam  theory,”  International  Journal  of  Solids  and  Structures,  Vol.  37, 
2000,  pp.  535-558. 

33.  Yu,  W.,  Volvoi,  D.,  Hodges,  D.  H.,  and  Hong,  X.,  “Validation  of  the  Varia¬ 
tional  Asymptotic  Beam  Sectional  (VABS)  Analysis,”  AIAA  Journal,  Vol.  40,  2002, 
pp.  2105-2112. 

34.  Yu,  W.,  Volovoi,  V.  V.,  Hodges,  D.  H.,  and  Hong,  X.,  Manual  of  VABS ,  2008. 

35.  MD  Nastran  R3  Release  Guide,  2008. 

36.  Kuhn,  P.,  Stress  in  Aircraft  and  Shell  Structures,  McGraw-Hill  Book  Company,  Inc., 
New  York,  NY,  1956. 


67 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OM&  No.  0704-01 8& 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing 
data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or 
any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate 
for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that 


notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently 
valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

18-06-2009  Master’s  Thesis 

3.  DATES  COVERED  (From  —  To) 

Jan  2008  -  Jun  2009 

4.  TITLE  AND  SUBTITLE 

Structural  Optimization  of  Joined-Wing  Beam  Model  with  Bend- 
Twist  Coupling  Using  Equivalent  Static  Loads 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Green,  Nicholas  S.,  LCDR,  USN 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Institute  of  Technology 

Graduate  School  of  Engineering  and  Management  (AFIT/ENY) 
2950  Hobson  Way 

WPAFB  OH  45433-7765 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

AFIT/GAE/ENY/09-J01 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory  Air  Vehicles  Directorate 

Attn:  Dr.  Maxwell  Blair 

210  8th  St  Bldg  146 

Wright-Patterson  Air  Force  Base,  OH  45433 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 
AFRL/RBSD 

11.  SPONSOR/MONITOR’S  REPORT 

NUMBER(S) 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


13.  SUPPLEMENTARY  NOTES 

This  material  is  declared  a  work  of  the  U.S.  Government  and  is  not  subject  to  copyright  protection  in  the  United  States. 

14.  ABSTRACT 

This  study  is  based  on  the  merger  of  two  separate  theories  to  further  the  efficiency  with  which  joined-wing 
structural  models  are  designed.  The  first  theory  is  Geometrically  Exact  Beam  Theory  (GEBT).  GEBT  is  a  small 
strain  beam  theory  which  is  capable  of  accurately  capturing  the  geometric  bend-twist  coupling  in  beam  elements 
that  are  experiencing  large  global  deformations.  This  is  crucial  to  the  joined-wing  problem  as  it  is  geometrically 
nonlinear.  The  second  theory  concerns  Equivalent  Static  Loads  (ESL).  These  ESL  consist  of  a  load  vector  that 
produces  the  same  nodal  displacements  and  rotations  as  those  computed  from  a  pure  nonlinear  analysis.  The 
ESL  displacements  and  rotations  are  then  used  to  calculate  ESL  stresses.  By  merging  these  two  theories  into  a 
single  structural  optimization  effort,  computational  cost  is  reduced  by  orders  of  magnitude  when  compared  to 
purely  nonlinear  response  optimization  efforts.  It  is  shown  that  the  final  design  obtained  by  the  optimization  is  the 
same  for  both  types  of  analysis.  The  final  result  is  a  much  simpler  model  than  a  detailed  finite  element  model 
of  the  joined-wing  aircraft  that  can  be  optimized  without  significant  loss  in  fidelity  in  a  fraction  of  the  time  required 
for  a  single  nonlinear  response  optimization  cycle  using  finite  element  analysis. 


15.  SUBJECT  TERMS 

Joined-Wing,  GEBT,  Equivalent  Static  Loads,  Optimization 


16.  SECURITY  CLASSIFICATION 

OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 

OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Swenson,  Eric  D„  Lt  Col,  USAF  (ENY) 

a. 

REPORT 

u 

b. 

ABSTRACT 

u 

C.  THIS 

PAGE 

u 

UU 

80 

19b.  TELEPHONE  NUMBER  (Include  Area  Code) 

(937)785-3636  x7479 

e-mail:  Eric.Swenson@afit.edu 

Standard  Form  298  (Rev.  8-98) 
Prescribed  by  ANSI  Std.  Z39. 18 


